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Analytical and Numerical Studies of Three- 


Dimensional Trajectories to the Moon 


A. B. MICKELWAIT* anno R. C. BOOTON, JR.** 
Space Technology Laboratories, Inc. 


Introduction 


Se PURPOSE of this paper is to present theoretical 
calculations of lunar trajectories from the point of 
view of the guidance problem. Although there has been 
much work in the past few years on the general problem 
of lunar vehicle orbits, little work has been presented 
that would be directly valuable in designing the guid- 
ance for a launch vehicle. The attempt here is to give 
results which will be immediately useful for present 
lunar flights. Since the vehicles assumed in this study 
are ballistic, the guidance problem is simply one of re- 
ducing the error at the end of powered flight, that is, at 
burnout. Errors at burnout are a function of the par- 
ticular design of a given vehicle, and consequently the 
analysis of the errors presented here is generalized to 
make it as useful as possible to all designs. Although a 
conventional two-dimensional approach is presented, 
the three-dimensional problems arising from a monopti- 
mum location of launch sites are also analyzed. To 
substantiate the conclusions of this paper, a large num- 
ber of machine calculations were made and are com- 
pared directly with the approximate solutions. 

Section (1) contains the conventional two-dimensional 
analysis with the usual decomposition of the actual N- 
body problem, the earth, moon, vehicle, sun, etc., into 
the more tractable succession of two-body problems. 
The analytic statement of the required accuracy to im- 
pact on the moon’s surface is then determined and com- 
pared with the actual accuracy needed. However, 
since the accuracy requirements given assume that the 
launching may be made directly into the plane of the 
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moon, this solution is not satisfactory for real launch 
sites which require that the vehicle be launched out of 
the plane or that the vehicle be turned into the plane of 
the moon, an extremely expensive alternative for pro- 
pulsion. Therefore, Section (2) examines the effect on 
the necessary guidance accuracy of launching out of the 
plane of the moon. This three-dimensional analysis in- 
dicates that the velocity accuracy requirements in- 
crease sharply as the plane of the trajectory is inclined 
with respect to the plane of the moon, and also that the 
accuracy requirements vary greatly for different launch 
times, both for given times of the day and given days 
of the month. For example, from middle latitudes the 
accuracy requirements coupled with propulsion require- 
ments increase so greatly during several weeks of the 
lunar month that it is impractical to attempt to come 
close to the moon. Methods for analytically determin- 
ing appropriate launch times and azimuths are also 
given. 

Section (3) is a brief application of the ideas pre- 
sented above to some of the special problems posed by 
attempting to establish a lunar satellite. 


Section (1) 


In the initial analytic approximation, the earth- 
moon-vehicle system will be treated as a succession of 
two-body problems. That is, the vehicle will be con- 
trolled solely by the earth’s field in the initial phase, 
and then pass into a terminal region where only the 
moon controls its motion. To further simplify the 
problem, it is assumed that the moon moves in a circu- 
lar orbit about an inertial earth; the influence of both 
the sun’s and the moon's mass on the earth’s motion is 
ignored. In this section it is also assumed that the 
vehicle at burnout is moving in the moon’s orbit plane. 
It can be placed there either by a fortuitous location of 
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Fic. 1-1. Planar geometry. 
Rw = distance from center of earth to center of moon 
w = angular velocity of moon 
7 = total time of flight 
r,@ = nonrotating polar coordinates 
a = burnout velocity 
89 = burnout flight path angle 
Ye = radius of earth 
hy = burnout altitude 
& = initial position of moon 


the launch site, or by expenditure of sufficient fuel to 
orient the vehicle’s velocity vector into the lunar plane. 
This maneuver is neither necessary nor efficient, but the 
analytic results are useful as a beginning when ap- 
proaching guidance problems in the nonplanar case. 

The notation used in this section is shown in Fig. 1-1. 
First, considering only the effects of the earth’s field, 
assuming it to be central, the vehicle’s equation of mo- 
tion may be written in an inertial frame of reference 
centered at the earth as follows: 


Sn it a Oe (1.1) 


9 


y? 
and r6 + 276 = 0 


where uy = GM x = 1.4076 X 10" ft.*/sec.?. 
ing 6 from Eq. (1.1) via the usual angular momentum 
integral obtained from Eq. (1.2), there results 


Eliminat- 


p= V2E + (Qus/r) — (J2/r’) (1.3) 
where /: is the vehicle’s total energy (per unit mass) and 
J is the vehicle’s angular momentum (per unit mass). 

Because we ignore the moon’s mass in the initial 
phase, the interesting class of orbits, which are true 
minimum energy trajectories because the mass of the 
moon is used to pull the vehicle to the moon, are also 
ignored. These orbits consist of multiple revolutions 
of earth before arrival in the moon’s vicinity. Mini- 
mum energy, as used here, means the minimum energy 
to reach the moon’s vicinity while moving in an ellipse 
and without considering the moon’s field. The differ- 
ence between the two definitions implies a change of 
burnout velocity of approximately 200 ft./sec. The 
class of orbits ignored pose guidance problems too dif- 
ficult for the primitive vehicle considered here so that 
this definition of minimum energy is justified in this 


context. 
We are interested in the time required to go from a 
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radius, ro, to a radius, 7;. This is obtained immediately 
from Eq. (1.3) as follows: 


; { dr im dr , 
= _ = (1.4 
odd ro V2E + (2ue/r) — (J?/r?) 


Since J = r°6, we can obtain the total in-plane angular 


change, 0, from 


a " J dr 
§= — dr = ; = 
rn WF rr rV2Er? + Quer — J? 


Eq. (1.4) and (1.5) may be integrated in closed form 
(see Appendix) or integrated numerically. Using an 
earth radius of approximately 2.09 10’ ft., a mean 
circular lunar orbit of radius 1.261 X 10° ft., the flight 
time and in-plane angle can be obtained as functions of 
burnout angle, and burnout velocity as in Figs. 1-2, 
1-3, 1-4, and 1-5. Unless specifically varied, the burn- 
out altitude, 4p, and burnout angle, 8), are assumed to 
be 200 miles and 70°, respectively. These values are 
reasonable for a multistage, continuous burning vehicle 
of the type assumed here. Of course, fo1 detailed design 
work the interaction of 8) and 4 with the powered stage 
and subsequent ballistic flight must also be considered. 
Also, we limit the range of burnout velocities under 
consideration between minimum energy and approxi- 
mately 36,000 ft./sec. since we assume that excess burn- 
out speed will be converted into useful payload unless 
compelling guidance restrictions dictate otherwise. 
Because the effect of the moon has been ignored, the 
flight times and in-plane angles calculated from Eqs. 
(1.4) and (1.5) should be only approximately valid. 
The approximation should become increasingly poor 
as a minimum energy is approached and the vehicle 
spends more time in the combined field of earth and 
This is borne out by the curves 


(1.5) 


moon before arrival. 
labeled ‘exact results’? which represent exact numerical 
integration of the complete three-body problem. 

From Fig. 1-2, we see that the error in the analytic 
prediction of in-plane angle, @,* decreases from ap- 
proximately 4° at minimum energy to about 1° at 
36,000 ft./sec. From Fig. 1-3, we see that the predic- 
tion in time of flight improves even more rapidly as the 
velocity increases from an error of approximately 5 X 10° 
sec. at minimum energy to approximately 8 X 10° sec. 
at 36,000 ft./sec. This analytic model, of course, 
breaks down completely in the limited range of veloci- 
ties above exact minimum energy where the orbits may 
make numerous rotations about the earth before impact- 
ing on the moon. On the other hand, as the burnout 
velocity approaches infinity, the vehicle is effectively 
being aimed at a point moving in the moon’s orbit and 
the present approximation becomes exact. 

It is evident that the analytic results are inadequate 
for detailed orbit calculations unless burnout velocities 
are several hundred feet per second above minimum. 

*The in-plane angle, 6, as defined here is the total angle 
measured at the center of the earth between the initial radius 


vector andr = Ry». 
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¥ = Vehicle Velocity Relative to the Moon 
mM 
=, Impact Parameter 
ay, = Radius of Moon 
d = Vehicle Distance from Moon upon 


Entry into Moon's Field 





Fic. 1-7. Effect of the moon’s gravity on impact conditions. 


However, the model may be made sufficiently accurate, 
for guidance error analyses and design purposes, by 
turning off the earth’s field when close to the moon 
and calculating the correction due to the moon’s in- 
fluence alone. 

We now proceed to introduce the effects of the moon. 
To simplify the calculation in the vicinity of the moon, 
we will now ignore the curvature of the moon’s orbit 
geometrically as well as dynamically, and consider the 
problem in a moon-fixed coordinate system (see Fig. 
1-6). The velocity components in this coordinate sys- 
tem are 


= v0 
(1.6) 
y=e 


where vy, is the moon’s orbital velocity ~3,370 ft./sec. 
Using energy and angular momentum relations to com- 
pute the inertial injection velocity components. 


p= V2E + (Que/r) — (J?/r®) 
(1.7) 


The magnitude and angle of vehicle velocity with re- 
spect to the moon are then easily calculated from 


Umm = V x2 +y = V (on — ro)? + 7 
(1.8) 
tan~! [vy — 76/7] 


YW = tan! (%/%) = 


Now, the distance the vehicle crosses the moon’s 
orbit ahead of the moon in the absence of the moon’s 
gravity field is given by Ry(@ — wl — 6). This dis- 
tance is directly related to the “impact parameter,’’* 
b (see Fig. 1-6), by the relation 


* The impact parameter is the perpendicular distance between 
the moon’s center and the asymptote of the incoming hyperbola. 
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b = (cos Y)Ry(O — wl — A) (1.9) 


Thus, this model takes the conditions at the moon's 
orbit calculated without considering the moon’s field 
as the injection conditions for the hyperbolic orbit 
about the moon. 

The impact parameter, } (from the Rutherford scat- 
tering theory), is convenient since it determines the dis- 
tance of closest approach by recourse to the usual as- 
sumptions of energy and angular 
momentum in the moon’s frame of The 
small rotation of the coordinate system in the time in- 
If we assume the 


conservation of 
reference. 


terval under consideration is ignored. 
vehicle enters the moon’s gravity field at a distance, d, 
and if the vehicle’s velocity relative to the moon, 7,,y,, 
has an angle, e, with respect to the center of the moon 
(see Fig. 1-7), then 


b= dan < (1.10) 


Combining Eq. (1.10) with the energy and momentum 
relations, the maximum impact parameter allowable for 
collision to result is 


Omar = Ay V1 + (2ua/Am0ma") — (26a dimy”) (1.11) 


where ay, = moon’s radius, and py = GMy = 1/S1.- 


375un. Since d ay for any reasonable choice of d, 


we can approximate bq, by 











Maximum impact parameter as function of vehicle 
velocity relative to the moon. 


Fic. 1-8. 


Fic. 1-9. Angular difference as a function of burnout velocity. 
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(Bo = 70°). 


/ 2 ‘ 
Bnar = Ay V 1 + (2454 /A uma”) (1.12) 


The maximum allowable impact parameter, Bnaz, is 
plotted in Fig. 1-8 as a function of v,, ys. 

We are now in a position to relate deviations in burn- 
out quantities to changes in (@ — w7’), thence to changes 
in 6, and finally to distance of closest approach. The 
integrals of 7 and 6, Eq. (1.4) and (1.5), can be dif- 
ferentiated and then integrated in a closed form to ob- 
tain the first derivatives of (@ — w7) with respect to 
burnout velocity and flight path angle (see Appendix). 
This procedure is inadequate where 0/0v[(@ — w/’)| 
goes through zero (see Fig. 1-9) since higher derivatives 
play the dominant role there. In this region the results 
were computed from tabular values of (@ — w/)) asa 
function of velocity. The maximum allowable 6(6 — 
wl’) is computed from Eqs. (1.9) and (1.11) and then 
converted to maximum allowable errors in burnout 
quantities through the procedure described above. It 
should be pointed out that there is an element of arbi- 
trariness in the definition of a standard orbit with re- 
spect to which the allowable dispersions of vp and {8 are 
calculated. The standard orbit used here is such that 
the velocity relative to the moon, v,,,,, is directed at the 
center of the moon at impact. Because (@ — w/’) is not 
a linear function of vp for changes of vp as large as con- 
sidered here, the allowable spread of velocities is asym- 
metric about the standard. This asymmetry can be 
utilized to increase the total velocity tolerance, if the 
standard w is centered in the total allowable velocity 
spread and 7,,,, on the standard orbit is not aimed at the 
center of the moon. This choice of a standard, of 
course, would place more stringent requirements on the 
control of other burnout variables. For instance, the 
total allowable variation in 8) would necessarily decrease 
since impact locations of trajectories with perturbed 
values of 8) do not exhibit the same asymmetry as the 
velocity perturbed trajectories do. 

The results plotted in Figs. 1-10 and 1-11 show the 
strong dependence of allowable errors on the initial 
velocity, v. A comparison has been made with an exact 
computer study of this two-dimensional problem and 
the results of this study are presented also in Figs. 1-10 
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and 1-11. We note that the discrepancy decreases with 
increasing velocity; in fact, the discrepancy in calcula- 
tion of total allowable (év)) decreases from 10 ft./sec. at 
minimum energy to ~0 ft./sec. at v = 36,000 ft./sec., 
while the discrepancy in total allowable (66)) decreases 
from 0.7° to ~0° over the same velocity range. 

The effect of 6) on the allowable velocity error is weak 
for very low or very high values of vp and is strongest for 
values in the region where the effects of change in v on 
6 and w7 approximately cancel. The allowable velocity 
error that is maximum for each value of §» is found to be 
a decreasing function of 8). For example, the maximum 
allowable error of approximately 240 ft./sec. for the 70° 
value 8) decreases to 210 ft./sec. for 8) equals 80° and 
to 175 ft./sec. for By) equals 90°. The increase in allow- 
able velocity errors obtained with smaller values of 6) 
is not obtained without paying a price, since the velocity 
v at which these increases 
rapidly as {> is decreased. 

The variations with respect to the burnout altitude 
have been ignored so far. Variations in Ay are roughly 
equivalent to variations in wv, the equivalence being 
established through the requirement of energy con- 
servation at burnout.* 


maximum values occur 


Section (2) 


When the vehicle flies in the moon plane, a two-di- 
mensional flight, the same geometrical relationship be- 
tween the burnout point and the moon’s position recurs 
once a day. On the other hand, when the burnout point 









































* For example: since 0v/Oh¢ = —(pp/voro?), for v9 36,000 
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Fic. 1-11. Allowable velocity error for impact on the moon 
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lies out of the moon’s plane, the same geometrical rela- 
tionship occurs only once every lunar month, and then 
only if we ignore, as we shall, such effects as the varia- 
tion of inclination of the moon’s plane with respect to 
the equatorial plane. This fact has two consequences: 
one is that a different set of burnout conditions must be 
chosen each day to hit the moon; the second is that the 
sensitivity of the guidance system to error changes 
each day. Thus, for a three-dimensional flight, one 
problem is to determine the optimum day of the month 
for firing; the second is to find the correct combination 
of burnout azimuth, velocity, and flight path angle. 

The notation used in this section is shown in Fig. 2-1, 
which is a schematic representation of the surfaces and 
angles centered at the earth. Before analyzing guidance 
sensitivity, we must set up a procedure to determine a 
standard trajectory. The general procedure that we 
will use is to pick a set of burnout variables, vo, Bo, ho, 
that are consistent with booster capabilities and from 
these determine the in-plane angle, 6. The in-plane 
angle uniquely determines the required azimuth and 
time of day of firing for each day of the month. At the 
same time, the in-plane angle and a particular day of the 
month determine the angle between the moon’s plane 
and the vehicle’s plane, «, which is the single quantity 
necessary to determine the degradation of two-dimen- 
sional guidance tolerance to three-dimensional guidance 
tolerance. Therefore, we can calculate the variation in 
guidance sensitivity from day to day throughout the 
lunar month to aid in finding the optimum day for 
firing. 

Referring again to Fig. 2-1, the following relationship 
can be deduced for the required burnout azimuth, a» and 
6, as functions of \,, the burnout equatorial longitude re- 
ferred to the lunar descending nose, and \,,, the moon’s 
position in the moon’s plane referred also to the lunar 


descending node 


§ = cos~! } cos (Ay — ¥) cos [¢p + sin~'(sin y sin 7)] — 
sin (Ay — y) sin [@o + sin—! (sin y sin7z)] X 


V1 — (sin? \,/sin? y)} (2.1) 
a = sin = = p= ?] (2.2) 
sin 6 
! (tan \,/cos 12) 


where y = tan 


In Fig. 2-2, we see 6 plotted as a function of \, and Ay, 
for latitude 30°N. The in-plane angle is essentially a 
linear function of \, and Ay over the range of angles 
shown there. These angles were chosen since, as it will 
turn out, they cover the regions of practical design in- 
terest. 

Combining Eqs. (2.1) and (2.2) we see that for a par- 
ticular burnout latitude, do, and for a particular in- 
clination of the moon’s plane, 7, the burnout azimuth, 


ao, is fixed once we fix \, and Ay. Of course, d, varies 


through 360° every 24 hours and essentially determines 
the time of day of firing. On the other hand, A, varies 
only 12° to 13° per day and within an uncertainty of 
about 6° or 7° determines the day of the month. \, and 
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Aw are already physically related but the rapid varia- 
tion of }» compared to \y makes it possible to ignore 
this relationship for the first approximation and treat 
the two as independent. Referring to Fig. 2-3, we see 
ay plotted for latitude ¢) = 30°N asa function of \, and 
Aw Over a range of angles that turn out to be of practical 
interest to a designer. Also plotted on Fig. 2-3 are 
curves of ap for a fixed 0, that is, for a given set of burn- 
out conditions. These latter curves are obtained from 
Eq. (2.2). 

The first general goal of this section is thus achieved, 
namely, the proper firing azimuth and time of day of 
firing for a given set of burnout quantities and on given 
days of the month. To make these results more con- 
crete, we have converted the in-plane angle, 0, into 
velocity, choosing By) = 70° and iy = 200 miles. The 
results of this substitution are in Fig. 2-4 which shows 
the variation of burnout azimuth throughout the month 
for various velocities. 

Several interesting results can be seen in Fig. 2-4. If 
the burnout velocity is less than v) = 36,200 ft./sec., the 
vehicle cannot be fired into the vicinity of the moon 
during a certain period of the month. This time period 
is symmetrically located about \,, = 270°, 
from the Northern Hemisphere, that is, the position of 
the moon furthest north (see Fig. 2-1). We notice that 
this restricted region becomes larger as the burnout 
velocity is reduced (8) and ¢» remaining fixed) because 
the in-plane angle increases and it becomes more dif- 
ficult to fit in the angular sector defined by @. Of 
course, if the burnout angle, 8), can be varied,* the 
angular sector can be changed and other launch times 
But for any given >, the general 
At burnout velocities just above 


when firing 


become available. 
symmetry holds. 
minimum energy, about 2 weeks of the lunar month are 
excluded as possible firing days. 

For any velocity, the orbit with the maximum al- 
lowable azimuth occurs when the moon at impact is 
furthest south, that is, when A,, = 90°. Since orbits 
with burned azimuths near 90° take maximum ad- 
vantage of the earth’s rotation, the class of orbits with 
burnout velocities less than ~36,100 ft./sec. which im- 
pact when the moon is halfway between descending and 
ascending nodes makes most efficient use of the booster 
capabilities. 

On the other hand, there are several reasons why this 
may not be the optimum orbit. First, the location of 
various downrange equipment may dictate that the 
vehicle flies in a certain band of azimuths. Fig. 2-4 
shows that for velocities less than 36,100 ft./sec., orbits 
with southerly azimuths are impossible, while for 
velocities above 37,200 ft./sec., northerly azimuths are 
impossible. In addition, as we shall show later, the 
orbit making optimum use of booster capabilities is not 
necessarily the orbit making best use of guidance 
capabilities. 

Although Fig. 2-4 refers to 45°N, the results for 
45°S can be attained by translating the curves 180° 

* Since a 8p of 70° is about optimum, any change will probably 


result in a payload penalty. 
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Fic. 2-1. Geometry of earth to moon trajectory 
= burnout azimuth oo = latitude of burnout point 
= burnout velocity #= in-plane angle 


= burnout flight path angle Ry = radius of moon's orbit 


= inclination of moon plane to equatorial plane 

= longitude of moon from its descending equatorial plane 

= angle between moon plane and vehicle plane 

= longitude of burnout point from moon’s descending equa- 
torial node 


3 o,, : 
Burnout latitude, @9= 30 N; in-plane angle, 6 (deg); burnout 
longitude with respect to moon's descending equatorial 
node, ;, (deg - equivalent to hours of the day) 
= 121 deg 
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Longitude of the Moon from its Descending Equatorial 
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2-3. Burnout azimuth as a function of time of day and day 
of month 





2-5. Maximum allowable latitudes of burnout for various 
times of the lunar month. 
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2-6. Inclination between the planes of the vehicle’s orbit 
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If B, = 70 deg and h = 200 mi, then 
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Fic. 2-7. Minimum inclination for various burnout latitudes and 
in-plane angles. 


along the Ay, axis and reading the azimuths as angles 
from the south. To illustrate this north-south symmetry 
further and show the variation of excluded regions as 
burnout latitude is varied, we note the following result: 
on a given day of the month, and for a given set of two- 
dimensional burnout conditions, the furthest north 
(south) a burnout point may be situated is for that orbit 
oriented due north (south). This result may be shown 
by finding implicitly, in Eq. (2.2), the maximum @y for 
a fixed 6 and \,, and variable \,,. 


Using the condition that a) = 0, Eq. (2.2) gives for 
this maximum latitude. 
Pimar = 7 — 8+ sin! (sin Ay, sin 7) (2.3) 


Values of ¢»,,,, from Eq. (2.5) are shown in Fig. 2-5. 
There we see that the restriction in days of the month 
vanishes for any velocity if burnout occurs at a latitude 
within about 30°N and 30°S. On the other hand, if 
burnout is at the North Pole, the minimum allowable 
velocity is 36,900 ft./sec. and firing can occur only on 
one day of the month. 

Two additional facts must be considered before prac- 
tical use may be made of Figs. 2-4 and 2-5. First, it is 
more useful to correlate the appropriate burnout 
azimuth and the maximum allowable burnout latitude 
with the launch day, not impact day. To do this, the 
moon’s angular rotation during the time of flight is 
subtracted from \j,, and one obtains an angular position 
of the moon that can be immediately converted to a 
particular launch day. The other fact that must be 
considered is the latitude change of the vehicle during 
powered flight, which itself may extend 10° to 20°. 
Curves of burnout azimuth for fixed burnout latitudes 
thus refer to launch points with variable latitudes. 
Since the latitude adjustment for powered flight is 
largest for azimuths closest to due north, these orbits 
with velocities below 36,200 ft./sec. will have large 
the neighborhood of the excluded 


adjustments in 
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regions. Since, all other things being equal, the correc- 
tion for powered flight reduces the required azimuth, 
the excluded period will enlarge. 

Another general criterion useful for optimizing an 
orbit is to reduce its sensitivity to burnout perturba- 
tions. With this criterion in mind, we turn to a con- 
sideration of the angle between the vehicle plane and 
the moon plane, w. 

The strong maximum in two-dimensional velocity 
tolerance (see Fig. 1-11) occurs because as the burnout 
velocity increases above minimum energy, the sensi- 
tivity of in-plane angle, @, to velocity increments de- 
creases until near 36,000 ft./sec. and equals the sensi- 
tivity of time of flight to velocity increments. Since 
the two angular position errors, that of the vehicle and 
the moon, are in opposite directions, the orbit becomes 
insensitive to velocity perturbations. However, if the 
orbit is inclined to the lunar plane, the two angular 
deviations no longer lie in the same plane and no first 
order cancellation can occur. Thus, the larger the in- 
clination between the two planes, the more the velocity 
flight is insensitive to errors in By) (see Fig. 1-3), so 
that only the vehicle’s angular position error, A@, is im- 
portant in this case, and the miss as a result of errors in 
By should be roughly independent of the angle between 
the planes. 

Referring to the geometrical relations in Fig. 2-1, we 
can write at once a relation between 7, \,,, and X, as 
follows: 

; sin \, sin [d) + sin! (sin y sin 7)] 

si 2 = 2 = (2.4) 

sin y sin 8 (As, Ans) 

in which 6 may be treated as a function », and dy 
through Eq. (2.1). The results of Eq. (2.4) have been 
plotted in Fig. 2-6 for burnout latitude 30°N. If @ is 
given a fixed value, then a relationship exists between 
d, and \y, and we arrive at the curves of u for constant 6, 
also shown in Fig. 2-6. 

The curves in Fig. 2-6 for fixed @ demonstrate a useful 
result that for any given set of two-dimensional burnout 
conditions, the minimum in w occurs on a day of the 
month such that the burnout longitude from the 
northern hemisphere is —90° with respect to the lunar 
equatorial descending node. For the southern hemi- 
sphere the minimum wv, of course, occurs when \, = 
+90°. 
(2.4) by varying u with respect to \,, and ), and holding 
6 fixed. 
tion of \,, and \, when @ is fixed can be obtained through 


These results may also be obtained from Eq. 


The required relationship between the varia- 


Eq. (2.1). Substituting 4, = —90° into Eq. (2.4), we 
obtain minimum uw for various 6, 
“ sin (dp mee 5S = ° 
SIN Unin = : Northern Hemisphere 
sin @ 
(2.5) 
; —sin (69 + 7) . : 
SIN Unin = : Southern Hemisphere 
sin 6 
The results of Eq. (2.5) are plotted in Fig. 2-7. Asin 


Fig. 2-5, we see the appearance of maximum allowable 
burnout latitudes when velocities are sufficient'y low 
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In fact, from Fig. 2-1 it follows that when the minimum 
inclination between the planes is 7/2 for a given in- 
plane angle, then the orbital azimuth is due north in 
the Northern Hemisphere or due south in the Southern 
Hemisphere. Therefore, when the minimum inclination 


| is 7/2, the orbit is starting from the maximum allowable 





latitude for that velocity. Fig. 2-7 indicates that the 
minimum uw is zero for 
do} < +72 


since once a month the burnout point lies in the lunar 
plane, so that Fig. 2-7 is complete as indicated. 

The angle between the planes, uw, is important for de- 
termining the change in guidance sensitivity. Guidance 
sensitivities in the more general case considered here 
may be found in a manner analogous to that used in 
Section (1). That is, we calculate the orbit to the 
moon in two parts: The first controlled only by the 
earth’s field, and the second controlled only by the 
moon's field. The results of the first calculation at r = 
Ry are taken as the asymptotic conditions for the 
selenocentric hyperbola. 

First, then, we determine the position of the vehicle 
relative to the moon at the time its distance from the 
earth is equal to the moon’s orbital distance. Suppose 
deviations in burnout variables have caused in-plane 
angle deviations, A@, and time of flight deviations, A7. 
Fig. 2-Sa gives the exact geometrical relationship for this 
situation. 
the moon’s orbit over the arc, w7, and the angular 
separation of the vehicle from the z axis when r = Ry. 
Referring to Fig. 2-Sb, which shows our approximation 
in detail, we see that P, the distance from the vehicle 
to the moon at the time 7 = Ry, is given by 


However, we ignore the small curvature of 


P= Ry|[(AO — AwT)? + 4A0AwT sin?(u/2)|”%? (2.6) 


The first term in Eq. (2.6), Ry(A0 — Aw7>), is the two- 
dimensional orbital miss, that is, the distance measured 
along the moon’s orbit between the vehicle at crossing 
and the moon. As in the two-dimensional case, the next 
step is to determine an impact parameter, }, and then 
from it, calculate a distance of closest approach assum- 
ing conservation of energy and momentum in the 
moon’s frame of reference. Letting v,,,, be the ve- 
hicle’s velocity relative to the moon, we can reduce the 
problem to that of determining the angle, n, (see Sketch 
|) since the impact parameter, 6, is 

6 = 


P sin n 


By considering the geometric relations in Fig. 2-Sb, we 
arrive at the following expression for n 
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Fic. 2-8(a). Geometry of earth-moon-vehicle system for a non- 


standard orbit 
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Fic. 2-8(b). Details of geometry in the vicinity of the moon 
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where ~ is the angle between the vehicle’s inertial 
velocity, v;, and the moon’s velocity, v,,, and a is the 
angle between P and v;. For most velocities of prac- 
tical interest, it turns out that 7 is close to 90° and the 
P itself may be used as the impact parameter. 

For example, if the vehicle burns out at 30° N lati- 
tude and —90° longitude from the moon’s descending 
node and with v = 36,000 ft./sec., 8) = 70° and hy = 
200 miles, we find the following results as the vehicle 
enters the moon’s field of influence. The vehicle is 
moving almost radially from the earth and = 85 
Its inertial velocity is almost perpendicular to the 
moon’s velocity and ¢ = 84.5°. Taking into account 
the other factors in Eq. (2.7), we find that » = 86 
In other words, the vehicle’s velocity relative to the 
moon is almost perpendicular to the distance between 
the vehicle and the displaced moon, and, therefore, 
this distance can be used directly as the impact parame- 
ter. For lower burnout velocities the vehicle’s orbit is 
less inclined to the moon’s orbit and 7 decreases and P 
becomes larger than the impact parameter. The 
guidance tolerance for impact calculated using P as im- 
pact parameter is therefore a lower bound on the more 
exact calculation which takes y into account and com- 
putes a smaller impact parameter. This is not critical 
since, as we shall see, with decreasing velocity the lower 
bound approaches the upper bound, that is, the two- 
dimensional velocity tolerance shown in Fig. 1-11. 
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Fic. 2-9. Maximum allowable velocity error for impact on 
the moon (burnout time chosen to minimize angle between 
vehicle plane and moon plane). 
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Fic. 2-10. Maximum allowable angular error for impact on 
the moon (burnout time chosen to minimize angle between vehicle 
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Fic. 2-11. Total allowable velocity error for impact on the 


moon (burnout latitude, 45° N). 


navcazssvvnncensnevenasacocvoencrsvonnnvineaetanvanncenneeannaenenegnny ii Me HTT jeuceveecnengegnergenanstnit11 


AUGUST, 


SCIENCES 1960 


In order to compare the two-dimensional and three- | 


dimensional velocity guidance tolerances, we select a 
burnout latitude of 30° and a time of day and month 
such that for each velocity considered, the orbit has the 
minimum inclination with the lunar plane, u»i,. Then 
for each velocity we use the machinery set up in Section 
(1) for the two-dimensional case to compute the maxi- 
mum allowable impact parameter, and the angular 


sensitivities 


Aéd = (00/dv) Av, Awl’ = w(O7/dv) Av 
This is sufficient, together with the value of 2 in, to 
compute the maximum allowable velocity error, which 
is shown in Fig. 2-9. Near the guidance minimum point, 
the term (A@ — Aw/7/’) vanishes so that 
P > 2Ry,AO sin (u/2) (2.8) 
and, therefore, the peak of the velocity tolerance curve 
is most sensitive to inclinations of the orbit plane. As 
the velocity is lowered, the first term in Eq. (2.6) be- 
comes more significant and thus the three-dimensional 
tolerance approaches the two-dimensional. 

Also, in Fig. 2-9 the results of exact computer calcu- 
lations for the same burnout conditions are shown. It 
it is evident from the 10 per cent discrepancy between 
the two curves that the analytic model is quite ade- 
quate for design purposes. 

Since in Section (1) it was shown that the time of 
flight is insensitive to By (see Fig. 1-4), the three-dimen- 
sional impact parameter calculated from Eq. (2.6) 
should reduce, in this case, to the two-dimensional. 
Fig. 2-10 shows a comparison of exact machine calcula- 
tions at 30°N latitude and minimum inclination, and 
indeed the exact curve is within 10 per cent of the two- 
dimensional exact curve. The discrepancy with the 
analytic curve in Fig. 2-10 at low velocities reflects 
both the inability of the crude underlying model to 
account for all the curvature in the orbit due to the 
moon’s field and the nonzero sensitivity of time of flight 
to changes in fp. 

Considering the velocity sensitivity, again it is ap- 
parent from Eq. (2.4) that, for a given 0, increasing the 
latitude of the burnout point increases “ and so reduces 
the total velocity tolerance for impact to result. 
Furthermore, the inclination of the two planes is also 
increased if the burnout point is not halfway between 
the ascending and descending lunar nodes. Both of 
these effects are demonstrated in Fig. 2-11 which gives 
the variation of velocity tolerance throughout the 
month for several burnout velocities. 

It can be noted from Fig. 2-11 that the velocity 
tolerance for the orbit with burnout velocity closest to 
guidance optimum, 36,000 ft./sec., has the largest per- 
centage decrease throughout the month while the lowest 
velocity is least affected. We also note that the three- 
dimensional velocity tolerance for lunar impact does 
not possess the symmetries with regard to \,, already 
noted for burnout azimuth. The cause of this is that », 
itself, takes one lunar month to complete a cycle. On 
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the other hand, similar curves will apply to 45°S if 
the \,, axis is translated 180°. 

We are now in a position to estimate what values a 
designer would be most likely to choose for burnout 
variables. If sufficient propulsion capability exists and 
velocity control is a problem, the total probability of 
a successful impact may be maximized by choosing a 
burnout velocity of 36,000 ft./sec. (see Fig. 2-11), the 
day of impact with Ay, = 34° and time of day with X, 
= —90°. From Fig. 2-4, we see that burnout azimuth 
for this velocity and time of launch is about 74°. Of 
course, if 8) other than 70° or iy = 200 miles were 
chosen, the numbers would change somewhat but the 
reasoning would be the same. 

However, if the designer wished to maximize the 
payload for a given booster configuration or if the 
guidance system had more difficulty controlling burnout 
flight path angle, 8, than w, then he would choose a 
burnout velocity of ~35,450 ft./sec. with impact day 
90°, burnout when A, > —90°, and a burn- 
The velocity tolerance would be 


when Ay = 
out azimuth of ~6S8°. 
only slightly less than maximum for this velocity but by 
afactor of 10 below the tolerance for 7) = 36,000 ft./sec. 


Section (3) 


The previous sections demonstrated that a relatively 
simple analytic model is capable of accurate prediction 
of guidance sensitivities for orbits designed to impact 
on the moon. In this section we show that the same 
ideas can be applied to guidance problems arising when 
the orbit is designed to satellite the moon instead. The 
allowable velocity and burnout angle tolerances for a 
successful satellite are not the same as those for impact 
on the moon and are functions of additional design fac- 
tors as well. 

Before discussing the guidance sensitivities, we note 
some of the well-known results concerning the re- 
stricted three-body problem that are useful as a back- 
ground for this discussion. The earth-moon-vehicle 
system of differential equations has a first integral of 
motion, the Jacobi integral. In a system of rectangular 
coordinates rotating with the earth and moon with 
origin at the barycenter of the system, z axis along the 


axis of rotation, and x axis towards the moon, the 
Jacobi integral can be written as 
J =v — w(x? + y?) — 2(Ueg + Uy) (38.1) 


where v is the vehicle's velocity in the rotating frame of 
reference, w is the angular rate of rotation of the moon, 
and Uy», Uy are the earth’s and moon's gravitational 
potentials at the vehicle’s position. For a given value 
of J the vehicle is restricted to a region of space defined 
by Eq. (3.1) and the condition that v > 0. The par- 
ticular Jacobi surface of interest to us is the largest one 
which is still closed about the moon since this provides 
a sufficient condition for a stable satellite around the 
moon. Setting J = J» for this surface, then if the 
vehicle is inside the surface with J = Jo, it will remain 
inside the surface indefinitely. We are ignoring, of 
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course, the small effect of the sun on the vehicle's 
motion. 

In Fig. 3-1 we see the projection of this surface on the 
x, y planes. The projection on the x, z plane is almost 
identical. This surface was computed for mean moon 


parameters, 


Ry = 1.26 X 10° ft., w = 2.66 X 10% 


rad./sec., and My/M, = 1/81.345. Superimposed on 
this surface are projections of terminal portions of 
orbits with burnout conditions similar to those of Sec- 
tion (2), namely, 8) = 70°, 4» = 200 miles, and various 
burnout velocities. These orbits burn out at 30°N lati- 
tude and at a time of month such that Ay = —90° (i.e., 
conditions for minimum angle between orbit plane and 
vehicle plane). 

The point of interest about these orbits is that the 
locus of points corresponding to a given time on an orbit 
lie roughly along a straight line, and are almost linearly 
distributed along this locus according to the variation 
in burnout velocity. From Fig. 3-1 we see that the 
locus is inclined to the earth-moon line by about 30° in 
the moon plane and turns out to be inclined by about 
12° in the plane perpendicular to the moon plane. The 
linear displacement along the locus is about 1.2 x 10° 
ft./ft./sec. in burnout velocity change. 

If we desire only that the vehicle ultimately be in a 
stable orbit about the moon, then locating the position 
of the standard orbit is relatively simple. As has been 
pointed out elsewhere, we must provide an additional 
velocity increment after entering the surface defined by 
J to convert the vehicle’s selenocentric hyperbola into 


Direction of 


motion of moon 








Loci of points 
according to 
fixed time 





Variation of orbits inside the critical Jacobi region as a 
result of burnout velocity increments. 


Fic. 3-1. 
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Plot of an unstable orbit about the moon (projection 
of orbit in moon’s plane). 


FIG. 3-2(a). 
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_ Fic. 3-2(b). Plot of an unstable orbit about the moon [pro- 
jected perpendicular to the moon plane, z axis along moon’s axis 
of rotation, same orbit as Fig. 3-2(a)]. 
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Fic. 3-3(a). Plot of an unstable orbit about the moon (projection 
of vehicle’s orbit in plane perpendicular to moon’s orbit). 
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an elliptical orbit. Ignoring for the moment the moon’s 
influence on the locus described above, the standard} 
trajectory should be such that the maximum segment] 
of the locus is included in the Jacobi surface at the time 
of adding the velocity increment. This assumes the 
velocity increment is added in the simplest possible 
way, that is, according to a predetermined time. 

With this criteria for the standard orbit, the velocity 
and angular tolerances at burnout are easily calculated, 
For the conditions presented in Fig. 3-1, the maximum 
line segment is about 3 X 10° ft. and therefore, the 
velocity tolerance is +125 ft./sec. The standard burn- 
out azimuth and flight path angle are easily calculated 
using the methods of Sections (1) and (2). The allow- 
able tolerances on these quantities are obtained by ob- 
serving that the locus can be shifted approximately 5 xX 
10! ft. laterally and 5 X 10’ ft. vertically before there is 
a significant reduction in length contained in the sur- 
face, Jo. Converting the lineal displacements into an 
angular displacement, A@, we use the methods of Sec- 
tion (1) and (2) to relate 00/08 to the maximum 
allowable Aé. The 06/08 to be used here is that ob- 
tained without considering the moon’s influence. The 
result for the set of burnout conditions in Fig. 3-1 is 
68, = +1.2°. 

The consequence of requiring that the vehicle satel- 
lite, instead of impacting, is twofold: first, the burnout 
velocity tolerance for a successful flight at burnout 
velocities just above minimum energy is increased until 
it is comparable to the maximum tolerance for impact 
(see Fig. 2-9), which occurs at 36,000 ft./sec.; second 
the burnout angular tolerance is held comparable to the 
tolerance for impacting the moon at minimum energy 
(see Fig. 2-10). As a consequence, the two optimum 
guidance sensitivities for impact have been achieved at 
a single velocity which is itself nearly the least possible 
to reach the vicinity of the moon. 

The angle between the planes, 1, will play a role here 
as well as for impact orbits since tilting the vehicle 
plane relative to the lunar plane has the effect of also 
tilting the locus with respect to the x axis. Therefore 
there will be a day-to-day variation in velocity guidance 
tolerance of the type found for impact orbits (see Fig. 
2-11). 

This is, of course, a grossly simplified picture of the 
guidance problems. The orbit locus is significantly dis- 
torted in the neighborhood of the moon (see Fig. 3-1) so 
that the velocities corresponding to orbits passing near 
the moon require special treatment. Another problem 
arises since some of the orbits, after the addition of the 
necessary velocity increment, will intersect the moon's 
surface instead of satelliting. On the other hand, if 
criteria of indefinite capture is relaxed, the tolerances 
can be enlarged. Figs. 3-2(a) and 3-2(b) and Figs. 3 
3(a) and 3-3(b) give examples of orbits with J > J 
which make many circuits of the moon before even- 
tually drifting back toward the earth. The analytic 
treatment of the behavior of such orbits is difficult and 
perhaps best conducted on a digital computer. 
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Appendix—Evaluation of Time and Angle 20x108 ft 
Integrals 
The time to go from a radius, 7, to a radius, 7;, can be r 16 
obtained explicitly by integration of Eq. (1.4) as 
: . r ae 
{= F (1) = F(ro) 
where the analytical form of F depends upon the sign r 8 
of the energy, #. If E is negative (vehicle velocity less 
than escape velocity), ¢ 4 
| m Moon 
Fir) = : | von Dit: sacs T T T T T T T T ’ 
2E T SHE aT 4 ‘ 4 -8 12 16 -20 24\.28x108 ft 
ME ; ( Me + 2Er )| a 
si! : 
V -2E V un? + 2ES? ‘ 
If E is zero (vehicle velocity equal to escape velocity), 
12 
F(r) = ((J? + usr)/3ug’)] V —J? + Quer Earth 
If E is positive (vehicle velocity greater than escape Fr -16 Rotating Coordinate System 
velocity), 
’ L 20 Total Time Elapsed 
Fir) = (1 2E)[ V [2Er? + Quer — J? — (ur/V 2E) X 7 3.0 x 10© sec 
At =2x 10? sec 
. . / + 9 
In (ug + 2Er + V2E V2Er* + Quer — J*)] b -24 RR Thr 
The angular change is obtained by integration of Eq. 2s 
(1.5) as se 
6 = G(n) — G(r) L . => 
where Gir) = sin ( Mer — i ) Fic. 3-3(b). Plot of an unstable orbit about the moon [pro- 
r — < apm jection of vehicle’s orbit into moon’s orbital plane, same orbit as 
rV uz? + 2ES* Fig. 3-3 (a)). 
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The Strain Analvsis of 
Solid Propellant Rocket Grains 


M. L. WILLIAMS* 
California Institute of Technology 


Summary 


There is often a wide divergence of interest between the missile 
engineer responsible for thrust performance on the one hand, and 
the motor tube manufacturer concerned with chemical composi- 
tion on the other. 
in structural integrity is, therefore, frequently hard put to com- 
promise both aims efficiently in arriving at a configuration, par- 
ticularly because of the physical properties of the grain ma- 
terial. Whereas conventional engineering analysis usually deals 
with linear relations between stress and strain, grain materials are 
potentially nonlinear and, generally speaking, more sensitive to 
time dependence, including the effect of temperature. The 
background relating to viscoelastic analysis, including the theo- 
retical work of Alfrey, Tsien, and Lee, is covered with special 
reference to grain-like configurations, 

The practical aspects of the elastic-viscoelastic analogy is 
brought out, and illustrated for steady pressure and temperature 
The use of previously developed pressure 


The propellant grain design analyst interested 


loading on star grains. 
and thermal stress concentration data for propellant analysis is 
described, and the importance of physical property data and 
material representation by mathematical models is discussed. 

Finally, the needs for some presently lacking physical property 
data and a satisfactory failure criterion for grain materials are 
stressed in order to instigate additional work in this area 


Introduction 


Mo" PRESENT missile propulsion systems are based 
upon either liquid or solid fuels. Whereas the 
coupling between the structural behavior of the fuel 
and the missile performance is relatively weak for 
liquids, the converse is true for solids. In the former, 
the main interaction is through loads which the liquid 
and its vaporizing pressurization imposes upon the 
pressure-vessel container. One of the most important 
is structural stability, where the geometry such as shell 
thickness to diameter ratio and fineness ratio are 
parameters. Another example is fluid motion, or 
sloshing, which affects dynamic stability and control 
analysis. Further, ground handling loads require de- 
sign consideration in the sense that overall performance 
within certain environmental conditions must be 
satisfied. 

On the other hand, and without attempting a com- 
parison of the relative merits of solid and liquid fuels 
for any given application, an obvious fact should be 
recognized. Whereas liquids can sustain only com- 
pressive hydrostatic stresses, solids will not only resist 
Presented at the Structural Design for Hypervelocity Vehicles 
Session, IAS National Summer Meeting, Los Angeles, Calif., 
June 16-19, 1959. Revised and received December 9, 1959. 


* Professor, Daniel Guggenheim Aeronautical Laboratory. 


hydrostatic or dilatational stress (including tension), 
but will also resist shear or distortion. 

Hence, from a conceptual standpoint, fuel has a spec- 
trum of solidification, passing from a liquid through a 
slush to a viscous mass, and thence to the limit of a 
brittle 
problems are eliminated and others introduced as a 
The most important single effect 


solid. As might be expected, some design 
function of solidity. 
is that a structural failure of the solid caused by pres- 
sure, thermal or environmental loads can induce cracks 
and thus precipitate uncontrolled burning and _ lead 
to a catastrophic, or at best inefficient, behavior with 
an obvious failure to meet performance requirements. 

The purpose of this paper is to review some types 
of structural analysis problems which arise in the de- 
sign of solid-propellant rocket grains, and describe the 
techniques available to handle them. 


Physical Characteristics 


Before delving into the details of grain analysis, it is 
appropriate to describe briefly the material with which 
we are dealing, its processing into the final configura- 
tion, certain physical and mechanical characteristics, 
and the environment to which it will be exposed.f ' 

The two general types of solid propellants are the 
double base, containing essentially nitrocellulose and 
nitroglycerin, and the composite wherein small par- 
ticles of oxidizers such as ammonium perchlorate are 
mixed in an elastomeric binder. 

Double-base grain compositions are thermoplastic, 
softening with increasing temperature, and are formed 
by two processes—with or without solvents. When 
solvents are used in conjunction with the nitropoly- 
mers, the initial mixture of usually nitrocellulose, nitro- 
glycerin and solvent roughly in proportions 2:1:2 1s 
either extruded or cast into the desired shape, after 
which the solvent is removed. In the solventless 
procedure, the two nitrous elements are mixed together 
and then heated to approximately 150°F. to obtain 
the proper consistency for extrusion after which they 
are cooled and ready for use. 


7 T. L. Smith has presented a very readable background ac- 
count for composite propellants in Elastomeric-Binder and Me- 
chanical-Property Requirements for Solid Propellants, Memo. No 
20-178, Jet Propulsion Laboratory (CIT), January 1959. It is 
a pleasure to acknowledge also several valuable discussions with 
him. A more general and elementary treatment may be found 
in Rocket Propellants, F. A. Warren, Reinhold Publishing Cor- 
poration, 1958. 
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Fic. 1. 


Illustration of an internal burning-star grain. 


Composite propellant grains are also usually pre- 
pared in two manners, extruded or cast. In the former, 
oxidizer is dispersed through a rubber gum stock by 
milling and then extruded through a die of the proper 
shape. It is then cured and the outer surface covered 
with an inert material to inhibit burning around the 
outer periphery. In such a form, it is loaded into a 
metal tube container and bonded or mechanically 
attached in place and referred to as a cartridge grain. 
An alternate method of forming, which is a very 
attractive method for large-sized grains where ex- 
trusion equipment becomes unfeasible, is a cast grain. 
Here the oxidizer particles and liquid-fuel elastomeric 
binder are mixed and poured into the metal casing, 
which serves as a mold. Usually the mold also con- 
tains a central star-contoured mandrel, and the vis- 
cous mass is poured between the mandrel and tube 
wall to produce, after cure, an internal perforated 
grain which is integrally case-bonded to the enclosing 
tube (Fig. 1). 

Current practice recently has been toward the use 
of cast grains, rather than the extruded, which in a 
sense is unfortunate because it has proved easier for 
the chemists to produce the desired mechanical proper- 
ties in the latter. On the other hand, ballistic consider- 
ations have required that cast grains be used in spite 
of the difficulty in producing chemical compositions 
which must have a relatively high ultimate elongation, 
even though a low-modulus and low-tensile strength 
(50 psi) are acceptable. 


Loading Conditions 

A grain is subjected to various loadings during its 
life. One of the obvious is the pressurization induced 
upon ignition. As the internal surface of the grain 
begins to burn, pressure in the grain builds up to 200- 
1,000 psi in the order of one tenth of a second or less. 
The load is transmitted through the propellant and 
substantially resisted by the case with the important 
parameter being E,h,/Eh, i.e., the relative stiffnesses 
of the case and the propellant. Depending upon the 
specific star design, which is selected to produce a de- 
sired pressure variation during the burning time, more 
or less severe stress concentrations are present at the 


star points. Itis in these areas that high elongation is 
necessary to prevent the material from cracking before 
it is burned away. 

Another critical condition occurs during the curing 
cycle. When the mass is in the viscous state and cur- 
ing begins, the grain shrinks but no strains are induced 
until the gel point is reached and initial polymerization 
creates the molecular structural lattice. Even at this 
point, internal strains are relieved to a large extent by 
viscous flow. Subsequently, however, as it begins to 
set up, the dilatational strains build up and the relative 
This 


is particularly true as the bond to the case becomes 


freedom of the molecules is further restricted. 
solid. Then tension and shear strains become increas- 
ingly important, depending upon the rate of curing, 
energy the relative thermal 
expansion of the case and propellant, and the geometry. 


chemical release rate, 

After curing, the grain is stored and here subjected 
to gravity stresses as well as environmental thermal 
stresses. Frequently, a grain will crack under these 
conditions of relatively low elongation after successfully 
withstanding high strains, but for a shorter time. 
Here a reciprocity effect between low strain rave—long 
time conditions, and high rate—short time is evident. 
However, whereas qualitatively the ultimate stress, 
g,, Will monotonically increase with increasing strain 
rate, the corresponding strain at this stress, ¢,,, will 
increase for a while and then decrease after some critical 
strain rate which depends on the temperature. This 
behavior is indicated in Fig. 2, taken from Smith,' for 
a pure rubber material. 

To repeat, this critical strain rate for maximum strain 
changes with temperature. Fortunately, it has been 
theoretically and experimentally established that it is 
possible to relate strain rate (or time to failure at 
constant strain rate) to temperature empirically,’ viz. 


log ar = log(tr t,) = —K,(T = T,) [Ky T iz cciad T,)| 
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where fr is the time to observe some phenomenon at the 
absolute temperature 7(°K.), ¢, is the time at a refer- 
ence temperature 7,, 7, the glass transition tempera- 
ture below which the material is essentially brittle, 
and the A; are approximately universal constants for 
many polymers. This relation, therefore, permits the 
master curve (Fig. 2) to be immediately extended to a 
wide range of temperature for most composite pro- 
This temperature range is from slightly below 


pellants. 
It is assumed, however, that the 


7, to above 100°F. 
extent of cure, percentage of cross-linking and oxidizer- 
fuel wetting are unchanged in this process of changing 
temperature. Thus, data recorded at one tempera- 
ture may be extrapolated to another temperature with- 
out additional experiments. In such a fashion, there- 
fore, Smith has tentatively postulated that similar 
tests for ultimate strain and stress at various tempera- 
tures can also be correlated for propellant composi- 
tions. Qualitatively, therefore, on the basis of this 
equation, one expects that an increase in strain rate at a 
constant temperature would produce on, say, the per 
cent strain at maximum stress, the same effect as a 
decreased temperature at constant strain rate. The 
principal results, quoting Smith,' are 


1. Values of the ultimate properties at a high 
temperature (approximately 160°F.) indicate the values 
that should be found at lower temperatures at very 
low rates. 

2. An increase in strain rate at high temperatures 
increases both the ultimate elongation and the tensile 
strength. 

3. At low temperatures, an increase in strain rate 
gives a decrease in ultimate elongation, although the 
tensile strength increases. 


Returning to the other main loading conditions, one 
also has those resulting from gravity or inertia body 
forces during firing. These act over a shorter time 
than the similar environmental storage conditions and, 
hence, lead to a different viscoelastic response of the 
propellant. The resulting stresses are usually shear 
critical on the propellant-to-case bond, although in some 
cases it is conceivable that acceleration forces during 
flight would act sufficiently long to permit axial slump- 
ing of the grain. Such deformation of the burning 
surface can be deleterious if it changes the effective 
port area of the motor nozzle. 

With this brief description serving as background 
material, let us consider now the type of stress analysis 
to be used. At the lower temperatures, below the 
brittle point or glassy modulus temperature, 7,, the 
grain may be considered fully elastic and analyzed as 
such on the usual basis of small deformations. On the 
other hand, as the temperature rises above 7%, the 
grain, depending upon the strain rate, tends to become 
increasingly viscoelastic and, hence, by implication, 
time-dependent on its behavior. 

Fortunately, there has been considerable theoretical 
and experimental work at least in linear viscoelasticity 
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and, within certain important restrictions, it is pos- 
sible to extend the techniques of elastic theory. 


Mathematical Background 


As emphasized earlier, there are two characteristic 
deformations in a solid, the first relating to extension or 
dilation, and the second to shear or distortion. In the 
one-dimensional form these are written, respectively, 

T= UY (1) 
o = Ee (2) 
where / and yw are the Young’s and shear modulus 
respectively. In the more general three-dimensiona 
form for small deformations of an elastic body one has 


oi; = AVSi; + Que; (3 ) 


where \ and yw are the Lamé moduli.* If now the ma- 
terial properties are time and stress or strain dependent, 
there would be a more general relation between stress 


and strain, say 


O;[o(xz, t)] = O>[e(x;, ¢)] (4) 
where QO; and QO, could be algebraic and differential 
operators. For example, in shear, O; = 1 and O: = u 
as in Eq. (1); simple tension O, = 1 and O. = Fas in 
Eq. (2). 


Another simple one-dimensional extension would be 
to consider the stress proportional to strain rate as well 
as strain itself in which case Eq. (4) would become 


a(t) = Ee(t) + nlde(t) /dt} (5) 


and for this special linear form O, = 1 and O, = n(d/dt) 
+ E. Continuing with this example, known as a 
Voigt representation, if a constant stress, oo, is applied 
to a tensile bar, Eq. (5) can be integrated to give 

e(t) = op/E[1 — er "| (6) 
which is seen to approach the static strain at large 
This behavior is known as strain re- 


times (Fig. 3). 
tardation (creep) and gives an idea of how a time- 
dependent stress-strain law is written in a simple case. 

Unfortunately, in a rocket grain the stresses are not 
one-dimensional, or necessarily linearly viscoelastic.t 
If, however, test data establishes the fact that a given 
material is linearly viscoelastic, then it is presently pos- 
sible to conduct, in principle, a three-dimensional anal- 
ysis based upon the work of Alfrey,* Tsien,® Lee,® and 
others. 

In his fundamental paper, Alfrey established that the 
stress distribution in an isotropic incompressible linearly 

. 7; fori Aj; uw E/2(1 + v), a 


€ = y;;/2 and o;; 
€1 + €22 + €33, and 6;; = 1 for 


vE/(1 + v)\(1 — 2v), 3d = «;; 
i = jand zero fori # 7. 

+ Suppose a time-dependent stress produces an associated time- 
dependent strain; then if doubling the magnitude of the stress 
holding the mode shape of the time variation the same also doubles 
the strain magnitude without changing the shape of its time 
dependence, the response is said to be linearly viscoelastic 











Vii 
pr 
we 


pli 
In 
res 


wh 
any 
rep 
ind 
pre 
pla 
A 
the 
wel 
the 
rest 
Lee 
Rea 
tech 
Alfr 
com 
Bi 
of tl 
tion. 
in t 
com] 
port 


and 


wher 


Then 


where 
coeffic 


Note 
bo, QO" 
LY. 
Also n 
is rela 





pos- 


ristic 
yn or 
1 the 


lus 
Ona 


has 


ma- 
lent, 
Tess 


(4) 


ntial 
= @ 
is in 


1 be 
well 


(5) 


dt) 
iS a 


lied 


(6) 


irge 

re- 
me- 
ase. 
not 
ic.t 
ven 
0S- 
nal- 
and 


the 
urly 


d 
for 


ime- 
ress 
bles 
ime 





SOLID PROPELLANT 


viscoelastic medium was identical to that in an incom- 
pressible elastic medium if the boundary conditions 
were prescribed upon the stresses alone. The signifi- 
cance of his contribution lies in the generality im- 
plied by the assumed linear viscoelastic relationship. 
In terms of Eq. (4), the stress and strain operators are, 


respectively, 


VU oO” 
QO; = An, (7a) 
, 2X, ot” 
N OQ” 
O. = b, (7b) 
é p> ot” 


which is the generalized version of Eq. (5). Thus, for 
any material whose time-dependent behavior could be 
represented as in Eq. (7) where a,, and b, are time- 
independent constants of the material, a method was 
presented for calculating viscoelastic stresses and dis- 
placements in many problems. 

After Tsien® extended Alfrey’s work by showing that 
the incompressibility restriction could be lifted, as 
well as including body forces, Lee® called attention to 
the fact that previous treatments implied the subtle 
restriction that Poisson’s ratio be time independent. 
Lee’s subsequent modification, which is in accord with 
Read’s’ general formulation although by a different 
technique, finally permits a complete generalization of 
Alfrey’s Theorem for an isotropic, nonhomogeneous, 
compressible media. 

Because of its relevance to the discussion, a statement 
of the appropriate relations follows, using Lee’s nota- 
tion. Suppose one writes the general stress strain law 
in two parts, using the deviator or distortion stress 
components, 5;;, i.e., the usual o;; with the hydrostatic 


portions subtracted off as 
Sij = oy (1 3) 6 ij OKx (S) 


and similarly for strain 


. 


Ci = €i3 — (1/3) bij;ex, (9) 


where the strains are related to the displacements by 
€i; = (1/2)[(Ou;/Ox;) + (Ou;/Ox,)] (10) 
Then the analogous relations for Eqs. (1) and (2) become 
P? [si] = Q* [eas] (11) 
R' [ou] = S* [ea] (12) 


where the linear differential operators with constant 
coefficients are defined as 


n 


p ra) 
p= a , ete. (13) 
2 bn 


Note again for example that in simple shear P? = 
bo, QO" = qo, so that si; = (go/Poles; = (Go/2h0)(2e:;) = 
uy. The ratio 90/20 is, therefore, the shear modulus. 
Also note that hydrostatic stress-strain elastic behavior 


is related through the bulk modulus k, such that R’ = 
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Y, 9° = So and thus o;; = (S0/ro)ei; = (3k)e.* Hence, 
k = $0/3rp. 

The principle used by Lee in establishing the visco- 
elastic analogy was to observe that if one proceeded to 
take the Laplace transform of Eqs. (11) and (12), there 
resulted relations of the form 


$ij(X,, P) = [O(p) P(p)| éij(xx, Pp) (lla) 


[S(p)/P(p)léi(xn, p) (12a) 
t t 


Gij(Xz, P) 


in which if the ratios of the transformed operators 
which are independent of coordinates x,—are identi- 
fied with “‘transformed”’ shear and bulk moduli, they 
can therefore be carried through the complete solution 
of the field equations for certain elastic stress boundary 
value problems. After solving for these transformed 
stresses by the usual theory of elasticity, the inverse 
Laplace transform is taken to produce the time- 
dependent stresses, strains, and displacements. 

Without elaborating at this point upon the details 
of the analysis, inasmuch as some example calculations 
for a pressurized hollow cylinder are given in the 
appendix, consider the immediate practical conclu- 
sion. Suppose we have an internal pressure, p;, acting 
on the wall or web of a long uncased hollow cylinder of a 
viscoelastic grain material. On account of the visco- 
elastic analogy, we can immediately state that because 
the boundary conditions are prescribed solely upon the 
stresses, the linearly viscoelastic stress distribution is 
identical with the elastic one. One way of interpreting 
this situation is to realize that the stress equations of 
equilibrium must be balanced at any time (viscoelastic 
acceleration is neglected). 

Specifically, if the inner and outer radii are a and }, 
respectively,’ 


o(r) = pit({1 — (b/r)*]/(AX? — 1)} (14) 
op(r) = pii{1 + (b/r)?|/(A? — 1)} (15) 


where \=b/a. The elastic displacement in this case is 


u(r) = [p; (1 + v)/E(? — 1)) X 
((1 — 2v)(r/b) + (b/r)] (16) 
and we observe that it depends upon the material con- 
stants; hence, if the material constants are time de- 
pendent, the viscoelastic displacement in contrast to 
the stresses will also vary with time. 
Suppose next, however, that this same cylinder is 
reinforced by a thin elastic case with material proper- 
ties /., v. and thickness h. Now the elastic solution 


becomes® 


: (b/r)? Pi — r*p’ 
or) = (p' — pi) = + — (14a) 
A* — |] A* — 

’ (b/r)? pi — r\*p’ , 

o(r) = —(p’ — pi) = = (15a) 
m= 3 \* — 
* Summing over the indices, o;; = oi + o22 + ¢ is) 3p, 
where p is the hydrostatic pressure and e;; = @; + @22 + ¢ 
3 = volume change. Hence, 3p = 3k(AV/V) or p = kR(AV/V 


+ These illustrative cases were prepared by L. D. Stimpson, 


based upon some collected formulas.® 
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b(1 + v) ee ee, 
u(r) = —— [C1 — 2»)(pi — A*p’)(r/b) + (pi — p’)(b/r) (16a) 
E(x? — 1) 
where using c as the outer radius of the case and \’ = c/a, 
: 2(1 — v)p, 
p = “ : / (17) 
i ts eee (2 — 1I)E(1 + v,) Dt (1 — Qe) 
= BIA aa aa ae : — Bea" 
(A? — WE + v) 
or for a thin case of thickness / 
: 2(1 — »v)p, - 
y= a (17a) 
‘ 5 (1 — v,*)Eb 
[1 + (1 — 2v)d?] + (2 — 1) 
(1 + v)EA 
and we observe that even in the elastic case the relative 
rigidities enter.* Hence, the stresses as well as the dis- series. It can be shown that the same behavior can be 


placements become time dependent. In this case, one 
must resort to the Laplace transform technique to 
calculate the viscoelastic response as has been done by 
Woodward and Radok." 

In summarizing the analysis method, therefore, it 
may be concluded that fairly direct means are avail- 
able, in principle, for calculating stresses, strains and 
displacements for a wide class of viscoelastic problems, 
providing that one is dealing with a linear material. 
Additional background material on the theory of visco- 
elastic analysis can be found in a recent survey paper 


by Lee."! 


Material Representation 


It is natural to turn next to the question of deter- 
mining the stress-strain equation for solid propellants. 
Classically, this is done by representing the material by 
models of spring and dashpot elements,* ® which are 
combined in order to exhibit the same overall behavior 
as the propellant. The simplest model is the linear 
spring, which is known as a Hookean and gives a direct 
proportionality between stress and strain. The spring 
constant is quickly recognized as being analogous to 
Young’s modulus. Another simple model is a single 
dashpot (Newtonian) wherein the stress is proportional 
to strain rate alone. A natural extension of these 
basic models is to either place a spring and dashpot in 
series so that strains are additive, or spring and dash- 
pot in parallel, in which case the stresses are additive. 
The former is known commonly as a Maxwell model, 
the latter as a Voigt model. These four models as 
well as a typical more general three-element model are 
illustrated in Fig. 3 along with the typical response 
curves. 

For many problems, it is expected that the Maxwell 
and Voigt models will provide a reasonable approxi- 
mation to actual propellant response for small strains 
over a limited time interval. However, when this is 
not the case, it is necessary to form models of either 
Maxwell elements in parallel or Voigt elements in 


* As a practical matter, for this particular configuration, it 
may be observed that if the case is strong enough to hold the 
internal pressure, usually Eb< E,h. 


represented by either form; however, the viscosity and 
spring constants appear in the stress-strain equation 
It is merely a matter of experi- 
Practical 


in different forms. 
mental convenience as to which one is used. 
consideration is that if there are more than three or 
four elements, computations using such models often 
become prohibitive except for very simple geometries 
which do not occur very often in grain analysis. 
However, a convenient way in which to generalize 
stress-strain response for a linear viscoelastic material, 
is to consider an infinite array of dashpots and springs 
either in Maxwell or Voigt form.'? For example, if 
we have an infinite array of parallel Maxwell elements, 
the overall strain of the model is the same as in each 
Maxwell element, and the overall stress is a summation 
By taking a limit in such a manner 
Maxwell elements becomes in- 


of the stresses. 
that the number of 
finite, while the total stress remains finite, an integral 
relation between stress and strain is generated. This 
allows us to write the model parameters as a continuous 
function, rather than as a combination of 
The experimental problem is 


discrete 
moduli and viscosities. 
then to determine this continuous function. 

Supposing, however, that one accepts (for discussion 
purposes) a two-element model; there is still the ques- 
tion of attaching numbers to 7, etc., or other constants 
a, and b,, in Eq. (7). Here, all the usual experimental 
problems occur as in the analogous simpler case of 
measuring Young’s modulus of metals as a function of 
temperature and evaluating any rate sensitivity. 

Conceptually, the problem is quite clear, and is some- 
what similar to the determination of the transfer func- 
tion of an airplane. Indeed, Blatz and Williams'* have 
explored the possibility of this type mathematical repre- 
sentation as an aid to the physical interpretation of ma- 


terial behavior. In terms of a transfer function for 


example, consider Eq. (lla). One could write for 
the a,,, 6, constant coefficient case 
é P(p) — , 
= = = KG(p), p = tWw (18) 


5 Q(p) 


and visualize applying a sinusoidal stress input of fre- 
quency w, measuring the strain response, dividing one 
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functions; viz., Q/P and S/R relating the distortion 
stress to distortion strain, and dilatation stress to 
dilatation strain. So far, the emphasis of most investi- 
gations has been restricted to a very simple special 
case—namely, uniaxial tension. Considering the pos- 
sible nonlinear behavior of the grain material, it re- 
mains to show that the transfer function remains 
the same 
under biaxial stress and hydrostatic pressure. 
is very little data at present for shear or torsion char- 


as would be predicted for linear material 
There 


acteristics because appropriate experiments and equip- 
ment have not yet been completely designed. 

While elastic behavior below the glassy modulus 
temperature is reasonably well known, material prop- 
erty determination for stress analysis purposes above 
7, is in an unsatisfactory state, particularly the deter- 
mination of a Poisson’s ratio, v. The task is made 
more difficult by the many load and environmental 
parameters to which the polymers are sensitive as well 
as the multiplicity of the materials now under active 
exploration for ballistic reasons. 


Failure Criteria 


Regardless of how well one may succeed in making a 
stress analysis, elastic or otherwise, a strength analysis 
cannot be completed until a failure criterion is estab- 
lished. 

Unfortunately, failure criteria for solid-propellant 
grains are in a very confused state, and have been too 
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Fic. 5. Contact pressure variation.!4 


confined to uniaxial testing. There are, first of all, 
conflicting requirements which dictate either or both 
stress or deformation criteria depending upon time. 
On the one hand, grain materials must withstand high 
short-time loads in contoured shapes but must not 
slump—like wet concrete—after long times. Briefly 
put, the usual failure criteria problems are encountered 


below 7’, and are compounded at the higher tem- 


peratures by being rate dependent. 
As a temporary expedient, one may invert the trans- 


form solution é;;(x;, ) to obtain 
€i; = €;:;(X,, C) 
and calculate the corresponding strain rate 
Rij = 0e,/Ol = KAx, 0) 


as a function of time. Then, knowing the strain and 
strain rate resulting from the applied loading, one could 
compare the ultimate strain data with test results 
analogous to the master curve given by Smith (Fig. 2) 
to determine a time-dependent strain margin of safety 
in the grain. The deficiencies such as arbitrarily 
selecting strain as against stress criteria and extending, 
without verification, uniaxial tensile failure to biaxial 
failure are obvious; improved procedures are required. 
In the elastic range, some progress has been made. 
It appears that a strain theory of failure may be ap- 
propriate, although whether it should be based on 
maximum normal strain or distortional strain is not 
clear at this time. To investigate the viscoelastic 
range, a series of tests should be carried out in con- 
junction with theoretical analysis for the biaxial stress 
state existing in a thick walled cylinder. This geometry 
is particularly favorable to analysis if the ends are re- 
strained and insulated so the heat flow is axially sym- 
metric and independent of cylinder length. There is 
only one equation of equilibrium connecting the two 
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unknown radial and tangential stresses. The second 
relation between these two quantities would, of course, 
be the failure criterion. Experiments might then be 
performed upon the cylinders at different temperatures 
and rates and various failure hypotheses for the viven 
material explored until one was found which areed 
to the desired Wise!* had started along 
these lines by exploring the normal strain hypothesis 


accuracy. 


for one material and finds agreement with test data in 
the low-temperature elastic range, but an investigation 
at the higher temperatures would be most desirable, 
particularly with a nonrigid case and unbonded ends. 
The object, of course, is to choose a geometry and 
loading condition which may be amenable to theoretical 
treatment and then run simultaneous experiments so 
that direct comparisons can be made. Once a tenta- 
tive hypothesis is established, more exotic geometries 
can tried. So far, at least, the hollow cylinder 
with its extensive background literature is the only 
type of experiment which seems to possess sufficient 
characteristic elements of the actual grain problem to 
be realistic. It is felt that in the current state, and 
particularly above the brittle point, the arbitrary 
prediction of failure in biaxial and triaxial fields of 


be 


strain or stress from uniaxial test data requires further 
substantiation. the 
hedral shear strain on an energy basis for triaxial failure 


Williams’ has proposed octa- 
and a cumulative damage law for generalizing the Smith 
failure curve to variable strain rates. 


Star Grain Analysis 


Having indulged in a good deal of descriptive material 
thought to be necessary in providing background for 
the analyst not particularly familiar with propellant 
grains or viscoelastic techniques, consider now a spe- 
cific illustration of what can be done for a fairly long, 





Geometry of a typical slotted grain 


Fic. 6. 
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grain as found by Ordahl and Williams.” 


Notwith- 
standing the pessimistic viewpoint taken in many as- 


internal-burning case-bonded configuration. 


pects of the work, it has proved possible to analyze 
certain shapes with a fair degree of precision, par- 
ticularly in the elastic range which is appropriate at 
low temperatures or high rates of loading. 

It is convenient to discuss three categories of loading; 


pressure, thermal, and environmental. 


Pressure Stress 

The first problems one naturally investigates for 
grain design are thick-walled cylinders under pressure 
loadings. This classical problem, for years a recurring 
one to gun barrel designers, is for most purposes ade- 
quately treated in standard texts* '® both with and 
without consideration of various end conditions of the 
cylinder. A collection directed specifically to grain 
analysts has already been mentioned.‘ 

Recently, Wise! collected the 
formulas for the completely restrained grain and pre- 
sented charts showing stresses not only as a function of 


has appropriate 


pressure, but also as they depend upon a steady-state 
constant-temperature differential across the wall. Of 
particular interest are the variations of hoop strain, 
€,, and pressure ratio p’/p; between the grain and case. 
These are shown in Figs. 4 and 5 and illustrate the 
strong effect of Poisson’s ratio, v. Considering the 
fact that most propellants are nearly incompressible 
rubbery materials certainly at the higher temperatures, 
and consequently have v = 1/2, it is observed that the 
hoop (and other) strains become vanishingly small. 
If the material is strain sensitive, this fact, as relates 
to failure, is quite important and would indicate that it 
was desirable to adjust the relative rigidities of the case 
and bonded grain such that the propellant was sub- 
jected solely to hydrostatic compression, and hence 
unstrained. Unfortunately, pressure loading is not 
always the overriding design condition and, hence, the 
aforementioned desirable situation often must be com- 
promised. The importance of the specific value for 
Poisson's ratio and its temperature and rate dependence, 
however, is evident. 

Turning to the star configurations, a useful applica- 
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grain by M. L. Williams.}9 


tion of the thick walled formulas arises from the fact 
that it furnishes convenient normalizing stresses to 
which the effects of changes of cross-section geometry, 
case bonding, or loading may be referred. This fact 
has been used by Ordahl and Williams” in studying 
internally slotted grains (Fig. 6). Design charts 
(e.g., Fig. 7) have been presented for the stress concen- 
tration factors due to internal pressure* which can be 
used in conjunction with configuration design charts 
for progressive, neutral, or regressive burning require- 
ments such as the type developed by Stone.'® 


Thermal Stress 


Still confining ourselves to elastic analysis but focus- 
ing attention on the thermal stresses, it is again found 
that the thick walled cylinder is the simplest con- 
A collection of appro- 
priate design formulas are referenced. On the other 
hand, and similar to the pressure loading case, when the 
wall thickness is nonuniform as in the case of an in- 
ternal burning star, analytical solutions are not readily 
obtained. Furthermore, numerical methods are not 
easily applied in areas of high stress concentration, 


figuration with which to work. 


and one therefore naturally turns to experimental 
methods. 

In order to complement the stress concentration 
factor design charts for star grains under pressure, 


* The suggestion for the initial experimental pressure setup is 
due to A. J. Durelli and his coworkers at Armour Institute, whose 
pilot experiments have been collected as Thiokol Chemical Corp. 
Report No. 6-56, March 1956, entitled Photoelastic Comparison 
of Stresses and Strains Associated with Different Perforation Shapes 
in Rocket Propellants. See also Society of Experimental Stress 
Analysis, Vol. II, No. 1, p. 55, 1953. 


1960 


Williams’® has obtained similar data for thermal stress 
(e.g., Fig. 8). The case considered is a free-ended 
grain subjected to a steady-state temperature difference 
across the web. 

Biot” and Muskhelishvili*! independently pointed out 
that certain classes of thermal-stress problems, to 
which the usual grain configurations belong, are anal- 
ogous to superimposing appropriate displacements on 
the boundary which are uniquely determinable by the 
prescribed steady-state temperature distribution or 
heat flux into the grain. In 1939, Weibel** exploited 
this analogy photoelastically and, after experimentally 
verifying the results by comparison with the analytical 
solution for the thick walled cylinder, evaluated the 
It re- 
quired therefore only minor modifications to Weibel’s 
experimental setup to test a rocket grain configuration. 
Tests for a series of slotted grains were completed 
wherein the maximum thermal stress at the star-point 
was found and then divided, for normalization, by the 
thermal stress which would exist in a thick-walled 
cylinder for the same temperature difference. The 


thermal stresses in a square thick walled tube. 


resulting thermal-stress concentration factors were then 
presented in charts as functions of the design param- 
eters.’® 

It is perhaps worth emphasizing at this point that the 
temperature distribution is usually a requirement for 
the subsequent thermal  stress.7 
Presuming the temperature or 
conditions are known at the surface of the solid, it is 


determination of 
heat-flux boundary 
not always an easy thing to determine the internal 
temperatures. Aside from analytical solutions which 
exist in certain cases either by direct quadrature or 
conformal mapping techniques for steady state, one 
again thinks of experimental determination or analogs. 
Thermocouples have frequently been used, particularly 
where transients are required, but for steady-state and 
boundary conditions a_ voltage 
This technique, using Teledeltos 
conducting fluid 
actually the one used in preparing to find the thermal 


relatively 
analogy works well. 
paper rather than 


simple 
medium, was 


stress concentration factors referred to earlier.’ 

An alternate semianalytical technique may also be 
used for determining approximate temperature dis- 
Minimum have been deter- 


tributions. principles 


7 Generally speaking, the three-dimensional classical thermal- 
stress boundary value problem requires a simultaneous solution 
of four coupled linear partial-differential equations in the three 
displacements and the temperature. It is only a consequence 
of the fact that thermal inputs do not, in most engineering prob- 
lems, accelerate elements of the elastic body sufficiently for those 
displacements to affect the temperature distribution. Thus, it is 
permissible to uncouple the four equations and solve the fourth 
(heat conduction ) equation independent of the displacements, by 
hypothesis, and, with the temperature field known, solve the 
remaining three equations for the displacements in three direc 
tions corresponding to the temperature. See the full equations 
in the Parkus?* paper. The only place in rocket grains where 
this separation might have to be further explored is under def- 
lagration—detonation conditions where the chemical energy 
terms would also have to be included in the foregoing set.?* * 
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mined for heat conduction by Parkus,** Biot,*® and 
Chambers:” that is to say, a functional representation 
of the temperature distribution may be assumed, say a 
power series with unknown coefficients which values 
are determined by setting the first variation of 
the heat conduction integral with respect to these un- 
knowns to zero and solving the resulting algebraic 
equations. For those acquainted with structural 
analysis, the method is identical in principle with the 
Rayleigh—Ritz procedure. 

Transient thermal stresses may also be important, 
although generally not during the burning because the 
rate of burning is faster than the rate of heat diffusion. 
On the other hand, during grain curing or environ- 
mental cycling, time dependent temperatures are needed 
for subsequent thermal stress calculation. Geckler* 
has charted certain temperature distributions for hollow 
cylinders and Nichols,** and Nichols and Presson” 
have numerically determined transient temperatures 
during curing cycles including the heat sources due 
to the chemical energy of the polymerization. 

Experimentally, it seems possible to measure transient 
thermal stresses in grain configurations by an extension 
of the photothermoelastic technique developed by 


Gerard and Gilbert.” 


Environmental Stresses 


The last major category of stresses to be mentioned 
are those associated with storage and flight perform- 
ance. One of the main advantages of the solid fuel is 
its state of continual readiness without extended pre- 
Hence, there is the implied storage 
which are 


flight check out. 
for possibly long times. For 
liable to creep under gravity body forces, a long- 


materials 
time analysis is necessary. Two cases are particularly 
pertinent insofar as analytical results are concerned, 
namely, the thick-walled cylinder lying on its side, or 
on its end. In both these cases, elastic analysis fur- 
nishes some guide lines, particularly when extended 
to the viscoelastic situation as illustrated in the appen- 
dix, if the material representation is known. 

There are, however, some important effects that 
need evaluation. The first of these relates to the 
finite length of the cylinder with accompanying end 
effects. When the grain is case bonded, the usual type of 
discontinuity stresses are introduced, usually in the form 
of localized shear. The practical solution to this com- 
plicated by the unwieldy nature of the elastic solution,” 
let alone any viscoelastic extension. Second, here is 
the possibility of thermal cycling during occasional 
stand-by check-out or shipping. A complex inter- 
action somewhat similar to the curing problem may be 
generated which also emphasizes the need for a failure 
criteria even if the viscoelastic analysis problem were 
solved. 

To some extent, of course, inertia loads are similar to 
gravity effects, although as indicated earlier the signif- 
icant time scale is usually much shorter. Here the 
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critical stresses are more sensitive to the detail design 
of the retention mechanism of the grain in the case, 
whether it be by bonding or mechanical means, par- 
ticularly at the grain ends where discontinuity stresses 
are present. 


Conclusions 


The main emphasis in this presentation has been to 
establish that there are certain direct methods pres- 
ently available for the structural analysis of solid- 
propellant grains. For this reason, the illustrations 
chosen have been simple ones so as not to obscure 
the general scope of the discussion. As such, they 
have been selected from series of cylindrical grains for 
which relatively simple solutions can be obtained 
multiplied by appropriate 
To ward off a complacent feel- 


analytically and stress 
concentration factors. 
ing, however, it should be obvious that designers are 
exploring many alternate configurations. The various 
problem areas described are certainly not all inclusive. 
For one example, the fineness ratio of the grain is just 
one parameter whose value makes a major difference 
in the practical calculation of a two- vs. a three- 
dimensional stress analysis. Another is the intriguing 
problem of the stresses at a continually burning and 
eroding boundary.*': ** Also, efficient internal ballistics 
usually require sharp corners in the grain to reduce the 
sliver fraction and, thus, increase the loading fraction. 
These are poor structurally. What is a good systems 
compromise? These are but a few of the associated 
detail problems. Vogel** has described some other 
possible grain configurations, and it takes little imagi 
nation to visualize the host of additional problems. 
The important point to remember is that elastic 
solutions are generally required before the viscoelastic 
problem can be treated. On this basis then, even 
disregarding such important restrictions as a knowl- 
edge of the material representation by mathematical 
models or potentially nonlinear viscoelastic effects, 
practical computational considerations will continue 
to rule out exact analysis of all problems of design 
interest. Important contributions must yet be made 
for approximate analysis techniques for viscoelastic 
media, no doubt utilizing energy techniques. ** 
Following the historical pattern of engineering analy- 
sis, the analyst must therefore compromise by con- 
tinually developing his intuition through the study of 
appropriate idealizations accompanied by judicious ex- 
perimental work. He can then permit himself edu- 
cated extrapolations to the specific problem at hand. 


Appendix 


A. Example of Viscoelastic Pressure Stresses 


Consider an infinitely long cylindrical hollow grain 
enclosed by an elastic case. The elastic solution [Eq. 
(14a), et seqg.| for an incompressible propellant (vy = 1/2) 
and a thin case under internal pressure in plane strain® 
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Consider the Voigt model representation with iso- 
tropic material constants of a viscoelastic behavior 


which evinces a time lag in shear 
(A-4) 


tiy = n(dyij/dt) + uyij (¢ = 7) 


In operational (Laplace) form with zero initial condi- 


tions, 
Ti = (MP + wy (A-5) 
As from incompressible elasticity theory 
Ty = BYy = (E/38) yy (A-6) 
Comparing Eqs. (A-5) and (A-6), we see the equivalent 
modulus becomes: 
E(p) > 3(np + pz) (A-7) 


This substitution for / is made into any of the initial 
equations to obtain the operational (viscoelastic) form 
for the respective stresses. Hence in Eq. (A-1) and 


(A-2) 
| + (>) | 4 4 hE, 
; cia eS 
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where upper sign applies for o,, and the lower for oy. 
Assuming a step input in pressure much like an ig- 
niter shock we have the appropriate transform 


pi = pi/p (A-9) 
and, hence, 
Grp) = RY (p + ao)/plp + (1/7)]f  — (A-10 
where 
k = —}[+(b/r)? — 1]/[(b/a)? — 1) fp; (A-11) 


hE, 
oor 


2(1 — »v,”)b 


The Laplace inversion of (A-10) gives the time-depend- 
ent physical stresses 


O(r9) = Rrj\ao — [ao — (1/r)Je”" }— (A-12) 


For example, when r = b, the pressure between the 


grain and case relaxes with a time constant 7, and gives 


a,(b) 1 — e WV? 
| »(a-) Ge) 
1 + 2(1 — »p,”) — 1 : 
a hE, 
The viscoelastic constant yz occurs in both the 


denominator and 7, while 7, the rate effect enters only 
in 7. Again it may be checked that when FE, > 
for a rigid case, r~! > o, and a,(b) ~ —p;. Note 
also the similarity to the elastic expression (17a) when 
it is specialized to an incompressible medium, i.e., 
y= 1/2. 

Finally, to obtain the suggested form for making the 
margin of safety calculation using Fig. 2, one computes 
for example ¢(a) and, by differentiation obtains the 


strain rate R(a) = (0/Ot)«e,(a) as follows: 


(= 4 1) hE, b? £4 
o9(a) a’ . 2(1 V-2)b itil a? u/s) a,(a) 
are (l—-e + Sistas - 
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b? hE, | b2 hE, 
Més(Q) a? 2(1 — v.*)ub _@,) ule,(a) oO =— 67) . go 
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= 1 
Pi ( b? ) hE, * b? 
“tate. se toa 
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These equations show that, after the initial step in pressure, the strain decreases from a peak value. 


, p; = ) . 
zo Be § 
a* 


The rapid 


pressure rise upon ignition has been presented in reference 8 where the rise was approximated by a ramp followed 
by a constant pressure and, alternately, by a damped sinusoid. In practice, this initial rise is nearly elastic. 


—p;, or hydrostatic compression, consequently 


* Note that for an infinitely rigid case (hE,/bE) ~ o and o, = of = a, 
e. = « = e, = O (Fig. 4 for vy = 


1/2). 








Si 


pr 


en 


de 





© Co. 
an ig- 


(A-9) | 


(A-10) 


(A-11) 








n the 
gives 


| 


1 the 
only 
> © 
Note 
when 
1.€., 


g the 
dutes 
; the 





apid 
wed 


ntly 





SOLID PROPELLANT ROCKET GRAINS 585 


B. Example of Viscoelastic Thermal Stresses 


Consider, again, the same cylindrical hollow grain 


with an elastic case. The elastic solution for the same 


uniform temperature rise of both the incompressible 


propellant and the thin case in plane strain is: 


o, = —p'} [1 — (a*/r*)]/[1 — (a?/b*)]}  (B-1) 
oy = —p'\ [1 + (a*/r*)]/[1 — (a?/b?)]}  (B-2) 
a, = —p’}1/[1 — (a?/b*)]} (B-3) 
where 
= | \ET | eee = | 
oo / 2(b? — a?) "ThE, 
(B-4) 


Consider this time a Maxwell model representation 





—6,(b) = p'(p) = 


Consider, as a simple temperature variation, a slow 
linear rise in temperature (somewhat compatible with 
the uniform temperature assumption) : 


T = k't; T(p) = (k’/p?) (B-10) 
Substitution into Eq. (B-9) gives 
b'(p) = k’/p[p + (1/7)} (B-11) 
where 
a (a — a,)uk (B-12) 
a~ + (1 2) bu 
= Ve . 
2(b? — a?) hE, 
a*u 
l 2(b? — a?) 


- = a? a . | 
— 93) — 
2(b? — a?) "Tae 1° 


The Laplace inversion of Eq. (B-11) gives again the 


simple lag function 
p(t) = kr(1 —e “”) (B-13) 


From the initial equations: 


o,(t) = —fkr{l — (a?/r*)|/[1 — (a? b*)] § x 
(4—e°“”) (B-14) 

o,(t) = —{kr[l + (a2/r?)|/[1 — (a?/b?)]} X 
(l—e “”) (B-15) 
o(t) = —}fkr/[l — (a*/b?)) {Cl —e “”) = (B-16) 


These equations show an exponentially increasing 
pressure from all sides as would be expected, since the 
Maxwell model has the capability of absorbing strain 
energy. 

The foregoing example was presented for purposes of 


demonstrating technique only. The assumption of a 


(a — a,)nup 


— p,?) 


| se 
2(b? — a?) 


of viscoelastic shear behavior: 


(dy;;/dt) = (1/u)(dr;;/dt) + (7;;/n) (1 #7) (B-5) 
In operational (Laplace) form: 
py = [(b/u) + (1/n) IF, (B-6) 
From incompressible elasticity theory: 
byy = (E/3)yvy = Ti (B-7) 


By comparing Eqs. (B-6) and (B-7), we see the equiva- 
lent modulus becomes in this case 


E(p) > 3nup/(np + wu) (B-8) 


This substitution for E(p) is made in Eq. (B-4) to 
obtain the operational (viscoelastic) form for the bond 
pressure. (The coefficients of linear expansion are 
assumed independent of these.) 


- T(p) (B-9) 
2 a~ 
- | np + - 


hE, 2(b? — a?) 


slowly rising uniform temperature is not very realistic 
for larger grains since the thermal conductivity of solid 
propellants is very low. In general, a transient tem- 
perature distribution exists which results in a diffusion 
rate for stresses and strains in addition to the time 


dependence resulting from viscoelasticity.’ 


C. Example of Viscoelastic Behavior Under Gravity 
Consider the hollow cylindrical grain supported by 
an elastically thick (rigid) case in an upright position. 
Neglecting end effects (or equivalently stated as an 
infinitely long cylinder) we have the grain subjected 
solely to one shear stress due to its own weight of: 


T(r) = —(p/2)|r — (a?/r)] (C-1) 


and a corresponding displacement (zero at case) of 


w(r) = — (p/2u){a? In (b/r) — [(b? — r?)/2|} (C-2) 


Introducing now a more sophisticated three-element 
viscoelastic model such that 
t((dr,-/dt)| + tre = m[t,(dy,2/dt) + Yr] (C-3) 
where 
T, = (nn/Um) = time constant for constant strain 


T, = Nm|(1/um) + (1/u-)] = time constant for constant 


strain 
Assuming an incompressible propellant (v = 1/2), the 
equivalent viscoelastic shear modulus becomes: 


m(Mte + Mm)P + Meltm : 
_ Inu Mm) ft (C-4) 


uA(7,p + 1) 
ii 
NmP + Um 


t.$ + 1 
Substitution of Eq. (C-4) into Eq. (C-2), and consider- 
ing a step input in gravity load such as perhaps being 


m 


suddenly uprighted 
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(C-5) 
Nn(Me 1 Mm)P TF bem 
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@(p) = k\(p + a)/plp + (1/7,)] (C-6) 
where 
p ¥ b b? — r’ 
k= -— a” In _ 
2(ue + Um) r 2 
Qo = (Um/Nm) = 1/7; (C-7) 


1, Tr = Mebm/ Nm He + fa) 


The Laplace inversion gives, for the axial slump, 


; s L L 2 »Xe 
wt) = a Jy In : -(*) }i- (4) |} x 
du, | r a b f 


t/tr 


1 + (¢m/Me) — (Um/Me)e 
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where it may be noted that, for long time, the elastic 
situation is approached and depends only upon the 
static material modulus p,. 
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Boundary-Laver Transition and Heat 
Transter in Shock Tubes’ 


R. A. HARTUNIAN,* A. L. RUSSO,** ann P. V. MARRONE*** 


Cornell Aeronautical Laboratory, Inc. 


Summary 


An experimental study is made of the wall boundary layer in a 
shock tube operated over a wide range of shock Mach numbers 
and pressure levels in air, including those for which real-gas effects 
exist. Transition distances are determined and correlated in 
terms of the transition Reynolds number based on a character 
istic length for this boundary layer. Data from independent 
shock-tube studies are also included in this correlation 

The results indicate a weak dependence of transition Reynolds 
number on shock strength up to moderate values of shock Mach 
number, followed by a larger stabilizing tendency. Comparison 
of these data with transition data obtained in the same manner 
in argon indicate that the increased cooling rates are largely 
responsible for the stabilization 

A dependence of transition Reynolds number on the unit 
Reynolds number is found at the lower shock strengths. Spe- 
cifically, higher transition Reynolds numbers are achieved at 
larger unit Reynolds numbers 

The phenomenon of transition reversal does not appear within 
the range of the experiments reported 

Laminar- and turbulent-flow heat-transfer rates to the walls 
of the shock tube are determined experimentally. The results 
of the heat-transfer measurements substantiate existing theories 
in both the laminar- and turbulent-flow regimes. 


Symbols 
= heat capacity 
C = ratio of product py at free-stream and wall conditions, 
Pelte/ Pwhtu 
( = specific heat at constant pressure 
h = enthalpy 
k = thermal conductivity 
Py = Prandtl number = cpyu/k 
g = heat-transfer rate 
Q. = defined in laminar heat-transfer relation as gwV x, see 
Eq. (2) 
Re = Reynolds number = p,(U, — U.)x/m 
St = Stanton number = q/p(U; — U.\(h, — h 
t = time 
7 = temperature 
U. = flow velocity in shock coordinates 
Ul’, = shock velocity relative to wall 
a = ratio of specific heats, c,/c 
\ = distance behind shock wave 
6 = boundary-layer thickness 
6* = boundary-layer displacement thickness 
7) = boundary-layer momentum thickness 
m = absolute viscosity 
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v = kinematic viscosity 

p = density 

T = shear stress 

¢ = ratio as defined by Eq. (3b 

Subscripts 

I refers to conditions ahead of shock wave 

é refers to region behind shock wave outside of boundary 
layer 

m refers to mean temperature for evaluating properties, 
see Eq. (3c) 

refers to recovery condition 
t refers to transition 
w refers to wall condition behind shock wave 


refers to properties of wall material evaluated at wall 
temperature behind shock wave 


Introduction 


hve EFFECTS of high heat-transfer rates on the 
transition of steady boundary layers from laminar 
to turbulent flow have been investigated experimentally 
by several research groups (for example, see references 
1 and 2). While the attainment of high heat-transfer 
rates is very difficult in wind-tunnel experiments, the 
boundary layer on the wall of a shock tube quite 
naturally experiences a high cooling rate. This occurs 
because the wall remains at essentially room tempera- 
ture during the test run, despite the extreme gas tem- 
peratures provided by the shock wave. In addition, 
the shock-tube wall boundary layer possesses many 
other characteristics that make it especially well suited 
for the study of many problems involving heat transfer 
and transition. Some of these characteristics are the 
inherent simplicity of shock-tube experiments, the ease 
with which significant parameters may be varied over a 
wide range, the fact that no models with their attendant 
complications are required and, in this connection, the 
absence of disturbances derived from the leading edge 
of models. Also, this boundary layer is amenable to 
theoretical treatment to the same extent that similar 
problems may be solved successfully for the flat plate 
in steady flow. Accordingly, an experimental study 
was conducted of the transition of the shock-tube wall 
boundary layer as well as the laminar and turbulent heat 
transfer rates to the wall over a wide range of shock 
Mach numbers and pressure levels. 

The boundary layer under consideration is that which 
develops from the foot of a shock wave traveling into 
a gas initially at rest. Such a boundary layer is 
sketched in Fig. l(a). Referred to coordinates fixed 
with respect to the shock [Fig. 1(a)], the shock-tube 
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Fic. l(a). Boundary-layer development in shock coordinates. 
Thermometer which moves with wall senses temperature history 
shown in Fig. 1(b) 





boundary layer is steady. From observation of the 
velocity profile, it is seen that this flow differs from the 
usual flat-plate case by the boundary condition re- 
quired at the wall. For shock-tube flow, the wall is 
in motion with the velocity of the shock wave. This in- 
troduces a nontrivial difference in the analysis of the 
flow. Viewed in the usual laboratory coordinates, 
the flow is, of course, unsteady. The characteristics of 
the boundary layer profiles are affected not only by 
the ratio of wall-to-free-stream temperatures and the 
Mach number of the flow relative to the wall, but also 
by the ratio of shock velocity to the velocity of the 
free-stream gas in shock-fixed coordinates (U,/U,). As 
the shock Mach number is changed in a given gas, 
all three of these parameters change. Therefore, the 
application of the present transition results to steady 
flows cannot be made directly. However, the develop- 
ment of a theory of stability for the shock-tube 
boundary layer would make it possible to determine 
the effects of significant parameters on the transition 
process which would have universal meaning. Regard- 
ing the experimental heat transfer results on the other 
hand, the existing theories make it possible to infer 
properties of high-temperature flows with high rates 
of heat transfer which are quite universal. For ex- 
ample, the range of validity of the assumption that 
pu is constant across the boundary layer may be 
assessed. 

Theoretical analyses for both the laminar and turbu- 
lent shock-tube wall boundary layer have been de- 
veloped and refined by several authors.*~® For the 
laminar boundary layer, the temperature measured 
at a fixed point on the surface of the wall is theoretically 
a simple step function, corresponding to passage of 
the shock wave. However, for turbulent flow, the 
wall temperature is expected to rise according to a 
fractional power of the time. Therefore, the tempera- 
ture history at a point on the wall indicates the time 
after arrival of the shock wave that the boundary layer 
undergoes transition. Figure 1(b) is an actual record 
of the wall temperature history. This record cor- 
responds to the flow situation of Fig. l(a). The resist- 
ance thermometer, being fixed to the wall, travels to 
the left of the figure at shock velocity in the coordinate 
system employed. Thus, as it crosses the shock wave, 
the gauge registers a step function in wall temperature 
[Fig. 1(b)]. The step function persists until the gauge 
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enters the region of transition to turbulent flow at which 


time the wall temperature increases with time. The 


plateau in the wall temperature record, which is indica 
tive of laminar flow in the boundary layer, is clearly evi 
dent as is the transition to turbulent flow. In addi 
tion, the heat transfer to the wall may be determined 
from the magnitude of the wall temperature rise. 
Therefore, from each record one obtains transition as 
well as laminar and turbulent heat-transfer data. 

Measurements of the wall temperature history behind 
a shock wave advancing in a shock tube were first re- 
ported at approximately the same time by Blackman’ 
and Chabai and Emrich. More recently, the results 
of many experimental studies, mostly at low shock Mach 
numbers, of the shock-tube wall boundary layer have 
appeared in the literature.°~ The data were obtained 
by various optical means and by measurement of the 
transient surface temperature. While most of these 
investigations were concerned primarily with heat 
transfer measurements, Chabai! and Asbridge!® also 
obtained transition data over an extensive range of 
conditions at relatively low shock Mach number. 

In the following sections, the present experiments, 
which extend the range of investigation to high shock 
Mach numbers, will be described, and the transition 
and heat transfer results will be discussed. 


Experiments 


The shock tube employed in this study was a con- 
ventional 1.5 X 2.5-in. rectangular cross-section tube. 
However, in view of the importance of surface smooth- 
ness on the processes of transition and heat transfer, 
a special 6-ft. parallel glass-wall test section was em- 
ployed. The glass-wall section consists of top and 
bottom walls of aluminum, lapped to a high surface 
smoothness, and glass sidewalls each of a single piece. 
The length was dictated by the requirement that all 
the air between the shock wave and transition originate 
within the section. In all but the highest shock Mach 
number runs, this condition was met. The glass also 
served to provide the proper backing for the surface- 
temperature gauges. 

Thin-film resistance thermometers were used to ob- 
tain the wall-temperature history from which the heat 
transfer and transition data were inferred. These 
gauges were also used to determine the wave speed. 
Resistance elements were located about 4 in. from the 
upstream and downstream ends of each plate glass 
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Fic. 1(b). Typical record of history of temperature at a 
point on the wall. Conditions of run were M, = 2.38; Pi = 


97.8 mm. Hg; Hg driver. 
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sidewall by deposition of a thin layer (about 1/10 micron 
thick) of platinum alloy and firing the elements and 
sidewalls together at about 1200°F. The calibration 
of the resistance thermometers and of the thermal prop- 
erties of the glass were determined by previously de- 
veloped techniques.’ The wave speed and _ wall- 
temperature histories obtained near both ends of the 
glass-walled section were reduced to provide heat trans- 
fer rates in the established manner.” Real-gas 
tables'*~*” were employed in evaluating all flow prop- 
erties. Sutherland’s formula was used to determine the 
viscosity coefficient. 

Tests were conducted throughout the shock Mach 
number range of approximately 1.5 to 10, in air, cor- 
responding to wall-to-free-stream temperature ratios 
of approximately 0.8 to 0.07. In each case the tem- 
perature of the air and shock-tube walls before the 
shock arrived was room temperature. The pressure of 
the shocked gas was varied to provide a range of Reyn- 
olds number per foot based upon free-stream prop- 
erties, (p,/ue)(U, — U,), of 0.5 & 108 to 3 & 108 where 
(U, — U,) is the velocity of the shocked gas relative 
to the wall. The allowable pressure behind the shock 
wave was limited by the strength of the glass sidewalls 
of the 6-ft. test section. Therefore, at the high shock 
Mach numbers, where the pressure ratio across the 
shock wave is large, only relatively low initial pressures 
of the test gas were used. Accordingly, at the high 
shock Mach numbers, tests were possible only at the 
lower Reynolds numbers per foot cited above. Both 
hydrogen and helium were used as driver gases to cover 
the full test range. The hydrogen driver, as usual, 
showed slightly more shock attenuation than helium;*! 
however, attenuation over the length of the 6-ft. test 
section did not exceed about 5% in most cases regard- 


less of driver gas. 


Results 


Transition 

In order to correlate the transition data for the 
boundary layer behind an advancing shock wave, it is 
necessary to define a Reynolds number that character- 
izes the flow. As pointed out in reference 6, the char- 
acteristic distance should be that which a particle in 
the free stream of the shocked gas travels from its 
initial position before it reaches the transition point, 
that is, (U, — U,)/U, Usts where t, is the transition time 
after arrival of the shock at the gauge. The character- 
istic velocity is taken to be that of the shocked gas with 
respect to the wall (U; — U,), and the kinematic vis- 
cosity is arbitrarily evaluated at free-stream conditions. 
Thus, the transition Reynolds number is defined as 


(U, —U,)* Ud; 
Me U, 


Re, = Pp, (1) 


The experiments were conducted keeping the Reynolds 
number per foot approximately constant (at different 
values for each series) and increasing the shock Mach 


number. As pointed out earlier, this causes the wall-to- 
free-stream temperature ratio, flow Mach number, and 
velocity ratio U,/U, to change. Figure 2 presents the 
transition Reynolds numbers obtained in the present 
and other experiments as a function of 7,,/7.. Shown 
also on the left-hand ordinate of this figure is the pa- 
rameter U,/U,. Itis seen that the data correlate fairly 
well at the higher cooling rates (7),/7), < 0.2), but for 
low cooling rates, that is, for low shock Mach numbers, 
the data show a definite dependence on the Reynolds 
number per foot. The apparent scatter in the data 
for a given shock Mach number (constant 7),/7), and 
U’,/U,) can be demonstrated to be due to a systematic 
dependence on the unit Reynolds number. Specifi- 
cally, higher transition Reynolds numbers are achieved 
for higher unit Reynolds numbers. This phenomenon 
will be discussed in more detail later. 

The range of values of the transition Reynolds num- 
ber goes from about 0.5 million to 50 million, the higher 
values being obtained at higher cooling rates and higher 
U,/U,. Not shown on this plot are four points ob- 
tained from reference 22 which go from 23 to 50 million, 
with higher values being obtained successively for lower 
7/7... The cross-sectional dimensions and surface 
finishes of the walls of the shock tubes employed in 
these various experiments are widely different. In 
addition, the detection of transition is accomplished 
by both optical observation and transient surface tem- 
perature measurement. Considering the great vari- 
ations in experimental conditions and techniques, the 
degree of correlation is seen to be good. Of the optical 
techniques used, only two are similar. Gooderum'* 
and Asbridge'’® obtained interferograms of the bound- 
ary layer on the bottom wall of a shock tube; Mark"! 
inferred transition by observing (schlieren photographs) 
the degree of interaction of a reflected shock wave with 
the boundary layer which that shock had produced 
before its reflection at the end wall. 

Smith et al.*? were able to observe the position of 
transition of the sidewall boundary layers from schlieren 
photographs. Figure 4 is a typical photograph ob- 
tained in the latter experiments. It is seen that a re- 
gion of smooth flow exists for some distance behind the 
shock wave, except for waves produced intentionally 
by inserting small plugs in both the top and bottom 
walls and accidentally by the inherent roughness on 
the steel walls of the shock tube. The beginning of the 
bubbly flow is interpreted as the onset of turbulent 
flow of the boundary layers on the glass sidewalls. 
That this turbulence is restricted to the sidewall 
boundary layers, and does not extend throughout the 
entire cross section of the shock tube, is deduced by 
the fact that the shock waves visible through the 
turbulent region are still quite regular. Also, by fol- 
lowing some of these waves, it may be seen that their 
angles change abruptly at some point. This is taken 
to be caused by the entrance of the shock wave into 
the sidewall boundary layers. The position of the 
interface between the driving and the driven gases is 
calculated to be a considerable distance behind the 
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Fic. 2. Boundary-layer transition correlation in terms of transi- 
tion Reynolds number. 


transition position. Some caution should be exercised 
in using the actual quantitative results of those ex- 
periments, in that considerable shock attenuation was 
reported to be present, and also, no special precautions 
were taken to achieve a long smooth surface. In addi- 
tion, large favorable pressure gradients were measured 
behind the shock wave, which might account for the 
greater stability of these data with respect to those 
obtained in the present experiments (see Figs. 2 and 3). 
Nevertheless, the qualitative consistency of the transi- 
tion data of those experiments with the correlation of 
the other data is very good. 

It is sometimes desirable to use boundary-layer di- 
mensions at transition as the characteristic length in 
order to take into account the effects of the parameters 
of the flow. In Fig. 3 the transition Reynolds number 
based on the boundary-layer thickness, evaluated 
theoretically** using the velocity profiles as viewed in 
laboratory coordinates, is plotted as a function of 7,,/7,. 
The data show the same trends as Fig. 2, including the 
unit Reynolds number dependence at low shock 
strengths. In addition, the data in Fig. 3 indicates a 
weaker stabilizing effect of higher shock strengths. 
This occurs because Re; is equal to the square root of 
the Reynolds number in Eq. (1) times a decreasing 
function of shock strength. In a preliminary pres- 
entation of the transition data,*4 the transition 
Reynolds number based on the displacement thickness 
using the velocity profile in shock-fixed coordinates 
was presented. This correlation showed a_ more 
decisive stabilizing effect of higher shock strengths 
than Fig. 2. Still another characteristic boundary- 
layer dimension is the displacement thickness using 
the velocity profile with respect to laboratory coordi- 
nates.* The Reynolds number based on this length 
shows a destabilizing effect at low shock strengths 
followed by stabilization for stronger shocks. Calcu- 


* The authors are indebted to Drs. S. Ostrach and H. Mark 
for helpful discussions which prompted more detailed considera- 
tion of the characteristic boundary-layer length. 


1960 


lations show that up to moderate shock strengths this 
Reynolds number is proportional to the Reynolds 
number which uses distance behind the shock, (U’/,), 
as the characteristic transition length. Since it is 
considered that this length is not the significant one 
(compared to the distance used in Eq. (1)), it is felt 
that the second-mentioned displacement thickness 
Reynolds number also plays.a subordinate role in char- 
acterizing transition. Until these points can be settled 
conclusively, it has been elected to present only Ke;. 
The major feature common to all the plots is a very 
weak dependence of transition Reynolds number on 
shock strength followed by a stabilization of the flow 
at high shock Mach numbers. Since not only wall 
cooling but also flow Mach number and the ratio 
U,/U, change as shock Mach number is varied, the 
trend of the transition data may be due to a combina- 
tion of these effects. For example, the flow Mach 
number varied from 0.5 to about 3.0, while L’,/U’, took 
on values of about 1.7 to 10. The influence of each 
effect acting independently can be assessed only through 
additional experiments where the parameters are 
changed separately, or it may be inferred through the 
development of a theoretical analysis for the stability 
of the shock-tube wall boundary layer. As a pre- 
liminary step in this direction, a separate series of ex- 
periments was run in argon in which the ratio 7,,/7, 
was varied from about 0.37 to 0.15. For values lower 
than this, it was found that heat transfer to the gauge 
by radiation from the argon obscured the transition 
data. Both the flow Mach number and the ratio U’,+ 
U, are considerably smaller in argon than for runs in air 
at the same ratio of 7,,/7.. For example, at 7,,/7, = 
0.35, the flow Mach number and U’,/U, in argon are 
0.8 and 2.4, respectively, while the values in air are 
1.3 and 3.5. The transition Reynolds numbers in 
argon are about the same as those in air at the same 
7,,/7,. On the other hand, when the argon data are 
compared to the air data at the same values of U/l’, 





Boundary-layer transition correlation in terms of bound- 
ary-layer thickness. 
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the argon transition Reynolds numbers are higher by 
some 50%. Since the flow Mach numbers of air and 
argon are about the same in this comparison, it is con- 
cluded that the higher transition Reynolds numbers 
are due to wall cooling. Extrapolating this trend to 
the region 7),/7, < 0.2, it is felt that the large sta- 
bilizing effect found in these experiments is due chiefly 
to the increased cooling rate. This could be checked 
by continuing the experiments in series of gases with 
different ratios of specific heats, y. 

In attempting to make a direct comparison between 
the transition Reynolds number, Eq. (1), of the shock- 
tube wall boundary layer and that obtained in steady- 
flow experiments on flat plates or cones, it should be re- 
called that the velocity and temperature profiles are 
different in the» two cases. Accurate quantitative 
interpretation of the shock-tube data can be made 
only when the appropriate stability theory has been 
developed. On the other hand, qualitative behavior 
of the two boundary layers may be expected to be simi- 
lar. Thus, it is expected that conclusions as to the 
existence of the phenomenon of transition reversal for 
flow over flat plates under conditions of extreme cool- 
ing’ may be reached from the shock-tube data. This 
statement is made notwithstanding the practical differ- 
ence between the leading edges in the two experiments 
and any effect this might exert on the phenomenon of 
transition reversal. Within the limits of the experi- 
mental operating conditions of the present study, al- 
though it has been noted that cooling exerts only a 
relatively weak stabilizing effect until high shock 
strengths are provided, no reversal of the effect of cool- 
ing rate on transition is observed. It has been pointed 
out that this phenomenon has been observed on a 
smooth cone only at values of the Reynolds number 
per foot greater than eight million.*° Inasmuch as 
the present data were obtained for values less than 
this, no conclusion as to the existence of the effect may 
be reached. However, the data of Smith et al.** were 
obtained at values of the Reynolds number per foot up 
to 13 X 10°, and these data show no transition reversal. 
Certainly, experiments over a wide range of shock Mach 
numbers at large Reynolds numbers per foot are desir- 
able. 

In both the present experiments and those of refer- 
ence 15, it appears that the transition Reynolds number 
is proportional to the square root of the Reynolds 
number per foot at the lower cooling rates. For higher 
shock Mach numbers, this dependence on unit Reyn- 
olds number is smaller for the data of the present 
experiments (Fig. 2). However, it is seen that if the 
Reynolds number per foot is greatly increased, as in- 
deed it is in the experiments of reference 22, the data 
show greater stability. This may be due, in part, to 
the fact that at higher unit Reynolds number the bound- 
ary layer is thinner, so that any adverse pressure gra- 
dient induced by the boundary layer would be smaller. 


Concerning the observed dependence of the transi- 
tion Reynolds number on the unit Reynolds number 





Flow 
DIRECTION 


Fic. 4. Schlieren photograph of flow in shock tube showing 
transition region Shock Mach number = 8.0, initial pressure of 
air = 3.00 psi. 


(pressure level), this same effect has been found in 
many steady-flow experiments on flat plates and cones 
in supersonic flow*: * and on airfoils in subsonic flow.” 
This phenomenon is not completely understood. A 
portion of the effect in steady flow has been associated 
with leading-edge bluntness.*’ *” Also, it has been 
suggested that the observed effect may arise from a 
dependence of the frequency and energy levels of dis- 
turbances in laminar boundary layers on the unit 
Reynolds number.* <A review plot of the unit Reyn- 
olds number effect has been presented in Fig. 17 of 
reference 31. 

In consideration of the possibility of leading-edge 
effects being responsible for the unit Reynolds number 
effect in the shock tube, it should be pointed out that 
the ‘‘leading edge’’ consists of the region of interaction 
of the shock wave and the boundary layer that de- 
velops at the foot of the shock wave. This interaction 
produces disturbances essentially different from those 
due to bluntness of a mechanical leading edge, and has 
been shown to be limited to a very small zone of the 
order of the boundary-layer thickness.** Therefore, 
while a leading-edge effect may also occur in the shock 
tube, the character and magnitude would be consider- 
ably different from that found due to a blunt leading 
edge. Additional research is required before the causes 
of the dependence of transition Reynolds number on the 
unit Reynolds number can be determined. 


Laminar Heat Transfer 

The temperature records |Fig. 1(b)| may be analyzed 
to yield heat transfer data which may be compared 
with theoretical predictions.‘~* The experimental 
wall temperature rise is compared with the theory of 
Mirels® in Fig. 5. The wall surface temperature rise 
varies directly as the square root of the pressure behind 
the shock wave. For shock Mach numbers less than 
approximately 6.5, real gas effects are not noticeable, 
and the theories of references 4—6 give comparable 
results. At higher shock Mach numbers, the results 
of Mirels should be used inasmuch as he obtained an 
exact numerical solution of the equations for Prandtl 
numbers of unity and 0.72. In this theory, chemical 
equilibrium was assumed to prevail, and the real-gas 
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. The laminar heat transfer rate from the shocked air 
ie is compared with the theoretical prediction® in Fig. 6. 
The heat transfer rate is given as: 
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Fic. 5. Glass wall surface temperature rise for laminar boundary- 
layer behind shock wave. 


properties of the shocked air were obtained from refer- 
ences 18-20. Since numerical methods were employed, 
it was not necessary to assume the product py to be a 
constant across the boundary layer and equal to the wall 
value. § The difference that this makes in the theoretical 
predictions of wall temperature rise and laminar heat 
transfer rate may be seen in Figs. 5 and 6. The experi- 
mental points are seen to be in agreement with theory 
when the variation of the product of density and vis- 
cosity across the boundary layer is considered. The 
laminar temperature plateau [Fig. 1(b)] was clearly 
identifiable in each run and did not appear to be affected 
by any attenuation phenomena. It might be pointed 
out that the application of more refined techniques 
has resulted in considerable reduction of the experi- 
mental scatter. 

At very high shock Mach numbers (about 10), the 
surface temperature rise obtained from the resistance 
thermometer was only a small percentage of that pre- 
dicted theoretically. The degree of ionization in the 
extremely hot shocked air appears to be sufficient to 
electrically short the platinum gauge. Further experi- 
ments at these high shock Mach numbers in both air 
and argon substantiate this. The authors recently 
showed** that this undesirable effect could be elimi- 
nated by the application of extremely thin coatings 
of a good electrical insulator (silicon dioxide). The 
success of this technique was proved in a series of ex- 
periments in highly ionized argon. Coated gauges 
gave heat transfer results in agreement with theory 
while uncoated gauges indicated considerably smaller 


surface temperature rises. 


where (kpc)x refers to the properties of the wall mate- 
rial. As may be seen, the heat transfer rate varies 
inversely as_ the Reynolds 
Accordingly, the representative parameter 


square root of the 
number. 
St WV Re was calculated over a range of shock Mach 
numbers, where Si is the familiar Stanton number and 
the Reynolds number is defined as (p,/u,)(Us U.)x, 
not to be confused with Eq. (1). The experimental 
results agree with theory and indicate that this param- 
eter is essentially constant over the shock Mach num- 
ber range, and is independent of pressure. 


Turbulent Heat Transfer 


The problem of turbulent heat transfer in shock-tube 
flow has been considered analytically® by extending 
empirical steady-flow flat-plate integral relations to the 
case in which the wallis moving. This theory expresses 
the heat transfer rate to the wall as: 


(h, — hy)» . 
qw = : ‘ (3) 


(U,—- UF," 


where the skin friction is given 














Fic. 6. Laminar heat transfer in shock tube 
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= 0.046 2 (1 =) >|" x 
= (). ) D _ 
aU? 5” U.) 6 


and 


Lm cs Pm vid = 
g = (3b) 
Me Pe 


The mean temperature is obtained from a mean en- 


thalpy given by: 
lim = 0.5(hy + hh.) + 0.22(h, — fh.) (3c) 


Numerical values of the ratio 0/6 (momentum thickness 
to boundary-layer thickness) are tabulated in references 
6 and 23. Experimental turbulent heat transfer rates 
behind relatively weak shock waves have been re- 
ported.’ Agreement between theory and experiment 
was found at a shock Mach number of 1.5, but not at 
a shock Mach number of 2.5. An interferometric 
study!* of the turbulent boundary layer at shock Mach 
numbers of 1.39 and 4.55 reveals that the velocity pro- 
file obeys the one-seventh power law assumed by 
Mirels at the lower shock Mach number, and in the 
transition region at the high shock Mach number, but 
it takes on the shape of a one-fifth power law in the 
fully turbulent portion of the boundary layer in the 
latter case. On the other hand, similar experiments’ 
indicate that for supersonic flow the one-seventh power 
law is valid, while for subsonic flows this is not so 

An extended range of comparison is provided by de- 
termination of turbulent heat transfer rates from tem- 
perature records obtained in the present study. For 
the most part, only data obtained for shock Mach 
numbers less than 5.5 were useful, because of the early 
arrival of the interface. However, earlier experi- 
mental data obtained by two or the authors at higher 
Mach numbers with almost completely turbulent flow 
are also shown in Fig. 7 together with two low Mach 
number points from reference 10. The theoretical 
curve was obtained from Eq. (3). | A simplified formula 
for the ratio of momentum thickness, 6, to boundary- 
layer thickness, 6, suggested in a private communica- 
tion from Mirels,* was employed at the higher Mach 
numbers. This formula, 6/6 = 0.317[1 — (U,/U,)I, 
was used for shock Mach numbers greater than 5.0. 

From Fig. 7, it can be seen that the theory and ex- 
periments are in reasonable agreement throughout the 
range investigated. The theoretical curve shows a 
decay at the higher shock speeds, whereas the experi- 
mental data does not. This could be due to the fact 
that the high Mach number experimental data were 
obtained on a rolled steel surface with a roughness level 
much higher than that of the glass. 

It is interesting to note that incompressible turbu- 
lent heat transfer relations, corrected for compressi- 
bility by evaluating state and property variables at 

* The authors are indebted to Dr. Mirels for his cooperation 


in connection with this research. 








Fic. 7. Turbulent heat transfer in shock tube. 


the mean temperature [Eq. (3c)], agree very well with 
. a - ‘ $/ o- 

experiment. The resultant value of St VRe = 3.7 X 

10~* obtained from this calculation is constant and 


does not decay at higher shock speeds. 


Conclusions 


The shock-tube wall boundary layer provides a means 
for studying transition at high rates of heat transfer 
and over a wide range of significant parameters. The 
clarity of the transient surface temperature records 
makes it easy to determine transition as well as laminar 
and turbulent heat transfer data. Since the equations 
governing this boundary layer may be treated theo- 
retically, this will permit significant heat-transfer 
parameters to be quantitatively inferred from the ex- 
perimental results. 

The transition data of the present experiments, and 
those obtained in similar shock-tube experiments, 
have been correlated by plotting the transition Reyn- 
olds number based on the distance which a particle 
of the free stream travels with respect to the wall 
before it reaches the transition point as a function of 
the ratio of wall temperature to free-stream static 
temperature. The results indicate that up to moderate 
shock strengths, transition of the boundary layer is 
weakly dependent upon the flow parameters, whereas 
a marked stabilization is observed for high shock Mach 
Since the wall-to-free-stream temperature 
U’. could 


numbers. 
ratio, flow Mach number, and the ratio U, 
not be varied independently in the experiments, the 
individual effects of each of these parameters on transi- 
tion could not be fully assessed. However, com- 
parison of the air data with those in argon indicates 
that much of the stabilization is attributable to wall 
cooling. 

The phenomenon of transition reversal did not ap- 
pear within the range of the experiments reported. 

An effect of unit Reynolds number (or, more properly, 
pressure level) on the transition Reynolds number was 
found at the lower shock Mach numbers. In particu- 
lar, the transition Reynolds number appears to be 
proportional to the square root of the pressure. At 
the higher shock Mach numbers, this effect becomes 
considerably smaller. 

The experimental laminar heat transfer results agree 
with the predictions of theory when proper account is 
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taken of the variation of the product py across the 
boundary layer. For shock Mach numbers in excess 
of 9.0, the ionization of the gas was sufficient to cause 
erroneous surface temperature indications by electri- 
cally shorting the platinum thin-film thermometers. 
It has been shown that thin coatings of dielectric sub- 
stances applied to the gauges overcome this difficulty. 

The experimental results on the turbulent boundary 
layer show good agreement with the existing theory. 
This agreement with theory through a wide range of 
shock Mach numbers seems to validate the convenient 
theoretical approach of employing an extension of com- 
pressible steady-flow relations to the shock-tube prob- 
lem. Also, it was found that incompressible turbulent 
heat transfer relations agree very well with experiment 
over the shock Mach number range investigated. 
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The Analysis of Redundant Structures by the 
Use of High-Speed Digital Computers 


WALTER J. CRICHLOW* ann GERNOT W. HAGGENMACHER** 


Lockheed Aircraft Corporation 


Summary 


Large-scale redundant structure analyses are currently feasible 


by the use of modern high-speed digital computers. This capa 
bility opportunely meets the urgent need to solve complex struc- 
tural problems which otherwise would be hopelessly beyond the 
capacity of the hand desk computer. However, the difficulties 
have now shifted from tedious hand computations to the prob 
lems of adequately representing the structure by a model and 
of the peculiarities of irregular geometrical configurations. 

A wide scope of problem types can be handled by a generalized 
program approach. Matrix formulation is used for the organi- 
zation Of input data and for handling data transfer in the large 
complex of subroutines, including the formation of equilibrium 
and continuity the deflections 


Simultaneous treatment of thermal expansions and _ plasticity 


conditions to final loads and 
is included 

The use of minimum-size redundant systems is emphasized, 
departing from the philosophy of cutting members to provide a 
statically determinate structure. Improved numerical accuracy 
and problem size capacity is gained for a given computer 

Examples are discussed ranging from simple plane-load diffu 
sion problems to pressurized fuselage cutouts and complex wing 


fuselage-shell intersection-type problems 


Symbols 


I general element load (axial loads N, bending mo- 


ments 7, shear flows gq). Specifically, the load 
distribution satisfying continuity 

Lo any load distribution satisfying equilibrium with ex- 
ternal applied loads. Called zero system 

Li zero system in equilibrium with dummy loads (for 
deflection computation ) 

L; self-equilibrated internal load systems, the redundant 


load systems. , Called internal systems 


Xj redundant multipliers to the L; 
U strain energy 
dS general elementary flexibility, such as 


(ds/EA ) for axial loads N, 

(ds/EI) for moments M, 

(dxdy/Gt) for shear flow q 
fictitious element loads—thermal or plastic parameter 


Pg 


Presented at the Structures Session, IAS National Summer 
Meeting, Los Angeles, Calif., June 16-19, 1959. 

*Group Engineer, Structural Statics Group, California Divi- 
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€ general strain (general elementary deflection, axial 
or angular ) 
é external deflection 
general sequence numbering subscript 


n degree of redundancy 

E number of elements 

R number of restraints or supports 

C number of independent equilibrium equations 
d number of dependent equilibrium equations 
G gridpoint 

ABC Cartesian reference system 


Matrix Notations 


f matrix of redundant (internal) systems 

F matrix of zero systems 

7 matrix of fictitious loads 

F matrix of final loads 

Gi matrix of dummy-load distribution 

D matrix of element flexibilities 

V matrix of coefficients uj; of continuity equations 
V, matrix of constants uj of continuity equations 
Y matrix of redundant scalars X; 


B matrix of deflections 
F is the transpose of f 


Introduction 


Sp PRINCIPLES of the analysis of redundant struc- 
tures have been well known for many decades, 
although not always clearly understood and often 
unjustly treated as a stepchild by the practicing engi- 
neers because of the large volume of computations in- 
volved in its treatment. However, the organization 
of the main procedure into matrix operations by 
Langefors' and Denke* has made large-scale complex 
redundant analysis practicable for solution by means 
of the high-speed electronic computer. The complexity 
now possible and the problems posed by rapidly prog- 
ressing field of contemporary and future designs have 
brought to light complications which the simpler cases, 
dealt with formerly, did not reveal. These are the prob- 
lems of the model representation of the structure, de- 
tails of the geometry, and their influence on the formula- 
tion of equilibrium equations, the degree of redundancy, 
the effects of nonlinear behavior and thermal expansion 
of elements, to mention some of the major points. 
Some of these problems have been discussed indi- 
vidually in various publications, but it seems appro- 
priate to make a summarizing presentation of the po- 
tentialities and problems of a general program for re- 
dundant analysis, with examples demonstrating its 
diversified field of application. Special attention will 
be given to the equilibrium conditions and the degree 
of redundancy, as well as the redundant load systems. 





596 JOURNAL OF THE 


Discussion 


Most problems of internal load distribution in a 
structure are indeterminate by use of the laws of 
statics (equilibrium) alone. The additional equations 
required to determine load distributions are provided 
by the consistency of deformations within the struc- 
ture. Three basic means of analysis are possible: 

(1) Displacements are defined as unknowns and equi- 
librium of forces provides the equations for solution. 

(2) Forces are defined as unknowns and continuity of 
deformations provides the equations for solution. 

(3) A simultaneous formation of all equilibrium 
and continuity equations® is possible from which all 
the forces and displacements can be solved in one sys- 
tem. 

This last method has the intriguing quality of re- 
quiring no attention to the apparently strange question 
of redundancy. Its disadvantage is the much larger 
number of simultaneous equations for a given size 
problem than either of the other two methods. 

The method of forces [point (2)] is chosen for the 
programed analysis discussed in this paper. The 
theorem of Menabrea—Castigliano is generally the most 
versatile formulation of the continuity condition. Al- 
though this theorem is strictly valid only in the elastic 
region, the computational procedure derived can be 
modified for extension into the inelastic range. 


By defining the coefficients 
uy == J LL dS 

the system of continuity equations can be written 

EF) = ] U1 + U2X » + — 

J = 2 UnX1 + UeX2+.. 
UpX, + upXet.. 
The solutions X; can be determined by any convenient 
means for any set of constants “;), and the final loads 
obtained from Eq. (1). The deflections are obtained 
by the dummy-load method: 
flection, a corresponding dummy load is applied and an 
equilibrium system L,, computed (any type of zero-load 


distribution in equilibrium with the dummy). The 
deflections e due to the correct load distribution / are 


For each required de- 


then computed: 


¢ = JS Lp, Ld-S (6) 


(2) Nonlinearity and Thermal Expansion 


The theorem of the minimum of strain energy is only 
valid for strictly linear load-deflection characteristics 
of each of its elements. Nevertheless, it is possible to 
use the resulting formulation of the analysis (Eqs. 
1 and +) in the nonlinear range with certain modifi- 


cations. 
(a) The simplest possibility is to replace the elastic 
properties in dS by the secant modulus in the elements 


AEROSPACE 
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(1) The Basic Equations of Redundant Analysis 
The basic principles of redundant analysis are sum- 
marized in this text to obtain a common denominator 
for further discussions. 
The correct load distribution in an 7 times redundant 
structure is given by 


L = Lo + XL + XoL» + es + clin = 


n 


Lb+>DXL; (1) 
1 
Let the total strain energy in the structure be 
U = Sf (L?/2)dS (2) 


The integral is extended over the whole structure. 
The continuity conditions, in the form of the theorem 
of minimum elastic strain energy, require 
ou 
OX; 


9Q 


0 ee re re (3) 


This leads to the system of continuity equations 


JS L;-LdS = 0 (4a) 
= SLL) + >> X,L)dS = 
l 


S LjLdS + > X; SL;LdS = 0 (4b) 
l 


GJ=l...4...n 


and uo = SL,LdS 


at ayyXi... + unXn + Un = 0 
Hb ui X,... + UenXn + Un = 0 (5) 
H+ uy Xi... + UinXn + uj = 0 

which exceed the elastic limit: dS, = (ds/F;A). 


Since this secant modulus, leading to the correct end 
strain, cannot be exactly predetermined, this procedure 
leads to successive approximations starting with an 
elastic analysis. After each change of the flexibilities 
dS, the terms uj; and uj in Eq. (5) must be recomputed 
and the reduction of the system of continuity Eqs. (4) 
or (5) repeated. 

Much less computing is involved in the following 
way: 

The term L-dS = (Ly 


deliver the correct deformation also by the introduction 


LX ,L,)dS can be made to 
of a fictitious load term P, retaining the original flexi- 
g £ 


bility dSp. 
two approaches (Fig. 1): e, = (L/E,A) (secant 
modulus), or e, = [(L + P)/E,A] fictitious load P 


The continuity Eq. (4) is 


Comparing the strain ez, expressed in these 


and elastic modulus £p. 
modified into Eq. (7). 


S LAL + P)dSy = 0 (7) 
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The fictitious load ? must be omitted from Eq. (1) (the 
equilibrium equation). In this approach, the coeffi- 
cients u,;; of the continuity Eq. (4) and its forward re- 
duction remain unchanged; it is only necessary to com- 
pute the solutions for a new set of constant terms “jo = 
J'L;-P-dSo with each iteration. This reduces the 
numerical work considerably. 

The thermal expansion can be treated in the same 
way by a fictitious load 7) = e7-/y)-A in each element 
involved, except that this involves no iteration. 

(b) A more sophisticated type of approach can be 
formulated using a combined load-strain diagram which 
includes the thermal strain e7 [shifting the load axis 
(Fig. 1)]. The procedure uses the tangent (#7) to 
the load-strain diagram in the flexibility dS;, plus a 
fictitious load 7p which is the intercept of the chosen 
tangent line with the shifted load axis (for combined 
yielding and thermal expansion).* The total strain 


e =e, t+ er = (L+ Tp)(1/ErA) (8) 
Thus, the continuity Eq. (4) is written 
S (L + Tp)L;-dSp = 
S (Lo + Tp) LjdSp + ox, JS L,LdSr = 0 (9) 


but Eq. (1) still reads 


£oh* SSK, 
1 


The deflections e must be computed from Eq. (10): 
e= JS Ly(L + Tp)dSr (10) 


The strain according to Eq. (8) [with L, 7p and Ep 
from previous solution of Eq. (9)| is used to find 7p 
and E, for the next iteration step. 

In this method again, both the coefficients u;; and 
constant i), and a complete reduction of the equation 
system (Eq. +) must be recomputed. A mathematical 
development of a similar approach is given by Denke.’ 

A simplified variation of this method is achieved if the 
slope Ey is not the tangent but an adequate approxima- 
tion to the stress-strain diagram in the plastic region, 
with only one change of flexibility (dS;) after the initial 


elastic analysis. 


(3) Matrix Formulation 


The main operations of the redundant analysis® in 
matrix notation are 
Eq. (1), F= hyo t+ f-X (11) 
the coefficients of Eq. (4), 
uj; = V=f'-D-f 


i 


and the constant terms 


uj = Vo = f'-D-(Fy + Ty) (12) 


(including fictitious loads). 


* Tp is negative on the positive load axis of Fig. 1. 


di | t TANGENT Ey 
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Fic. 1. Yielding and thermal expansion 


Partitioning the elements (and consequently the 
matrices) into m groups, subscript r, and excluding 
cross coupling (such as Poisson’s ratio) among any two 
or more blocks, reduces considerably the size of matrix 
multiplications 

m 
v= Dd SDes, 
For any two groups a and 6 where such cross coupling 
occurs an interaction flexibility matrix D,, is defined, 


Vo= Df’ DA(Fo + To), (13) 


and the following terms must be added into V and Vp: 


AV a» = fa Darfr + Is Das Te _ | 
fa’ Darfr + (fa’Darfr)’ (14) 


AV cap = fa’ Dar( Fo + To)s + fy’ D’ an( Fy + T)e | 
Finally, Eq. (4) is then 
VX +V,=0 and X = —V—!-V (15) 


and the final loads F are obtained from Eq. (11). 
The deflections e are given by the matrix B 


B = Gi)’ D(F + 7») (16) 


The explicit inverse V~' of the coefficient matrix is 
generally not required, only the solution of the simul- 
taneous equations; that is, the product (— V—!- V9). 
The forward reduction of the V matrix into the tri- 
angular form is done first,j retaining the reduction 
factors. Hereafter, any V» matrix is reduced with 
these factors, and the unknowns are obtained by the 
back solution. Since rarely as many different load 
cases (i.e., columns of constants in V») have to be 
analyzed as there are redundancies, this approach in- 
volves considerably fewer operations and is faster than 
the true inverse. 


Partitioned Analysis 

For large-size problems of redundant analysis, a 
partitioned solution can be of advantage in computing 
time and accuracy. The basis is the following: 

All operations leading to and including the formation 


7 Gaussian reduction. 
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4M oSe2 
Gi oP a 
-“G P 
2 2 4p, / " 
ca Pr Ge o ee % 
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ELEMENT ELEMENT ELEMENT 


Fic. 2. Element definitions. 


of the matrices /, /), 70, D are assumed completed for 
the total structure. Then the structure is partitioned 
into any number of m arbitrarily sized blocks (6) each 
of aredundancy mp. 

The redundancies of each block are those internal 
systems wholly contained in each block. Thus, the 
corresponding columns of internal systems are taken 
from the matrix f to form individual matrices /, for 
each block, such that 

Vy = fr’ Dh (17) 
for two different blocks a and 6, however, 
Veo = fa Df, = 0 (18) 

The redundancies connecting the different blocks 
are those internal systems overlapping between two 
or more blocks. All those ”, connecting systems sepa- 
rated from the original matrix f form the connecting 
matrix f,. The total redundancy 

m 
n=n.+ > Ny 
The blocks are now analyzed first individually for the 
external loads /), for 7», and also for the unit connecting 
redundancies, f,.. For each block 


Ves = fo Dip. Oi... sm) (19) 


external loads: unit redundant systems: 


Von = fo’ D( Fo + To) Va = hp DF. (20) 
Xo = + ne Vn 5 Vow X o = = V, i ™ Vo (2 l ) 
Fy = Fo t+ for Xo fn = fet fo-Xen (22) 


F, and f,, are loads and the connection systems, respec- 
tively, which satisfy continuity within each block. To 
satisfy the continuity of the block connections, the 
unknown scalars X, of the connection systems must 


be determined: 


V. = fa Dia = f' Dias (23) 

Voc = fo’D(Fy + To) = fD(Fy + T.) (24) 
Xe= —V,-"- Vo (25) 

F = Fy + fa:Xe (26) 


(4) Structural Model 


Many structures for which an analysis is contem- 
plated are too complicated and involve too many ele- 
ments for an analysis as is. A greater or lesser amount 


of simplification and lumping usually is necessary to 
arrive at a model structure which represents reasonably 
well the actual problem, compromising between the 
desire for detail and accuracy, and the amount of 
effort and work to be spent on the analysis. 


Elements 

A large variety of structural problems can be handled 
using three simple types of structural elements to build 
up the model. The elements are joined in the grid- 
points, which are located by space coordinates in any 
chosen Cartesian reference system. The three types 
of elements (Fig. 2) can be defined in the following way: 

(1) Bar or stringer element of constant axial force 
N between two gridpoints. 

(2) Beam or bending element. The bending moment 
M in a gridpoint of a beam, or frame-like structure, is 
the result of three confocal forces in equilibrium, acting 
in three consecutive points of two adjacent bars. These 
would be pin jointed at the common gridpoint without 
the specification of the beam element. For the equilib- 
rium conditions, the reaction forces at the gridpoint 
are expressed in terms of a unit moment. 

(3) Shear panel. A quadrilateral, shear-flow-carry- 
ing element with tangential forces along its side. For 
the equilibrium program, the shear forces are concen- 
trated in the four corner gridpoints (Fig. 2). 

Examples of model structures built up from these 
elements will be discussed later in the report. 


Equilibrium Equations 

Equilibrium equations are formulated in the direc- 
tions of the chosen reference system. The coefficients 
are the three-component reactions from unit element 
loads, computed from the coordinates of the gridpoints. 
(See Appendix. ) 


Degree of Redundancy 

The exact determination of the degree of redundancy 
of the model structure, or any part of it, is of consider- 
able importance for the successful analysis by means of 
redundant internal load systems. 

The redundancy occurs when the number of available 
equilibrium conditions is insufficient for unique deter- 
mination of element loads and support reactions. With 
C the number of independent equilibrium equations, 
FE the number of elements, and & the number of re- 
straints (or external supports), the redundancy is 


n=HE+R-C (30) 


If the model structure has G gridpoints, then six equa- 
tions per gridpoint exist; i.e., Co = 6G is the total num- 
ber of available equilibrium equations. It is always 
possible that a certain number (d) of dependent equi- 
librium equations exist. Thus, C = Cy) — dare the 
equations which actually can be used to determine 
loads. 
Hence, 


» 


n=E+R—(Q—-—d) =E+R-— (6G —d) (81) 
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Of the six equations per gridpoint, three are axial and 
three are moment equilibrium equations. In a struc- 
ture with no bending elements, the three moment equa- 
tions per gridpoint are dependent. It is possible to 
use the bending moment in a gridpoint as one unknown 
(as shown in element definitions) and express the three 
force reactions in terms of that moment. The moment 
equilibrium in gridpoints is then already satisfied and 
only three axial equations remain per gridpoint. Thus, 
Eq. (32) 

n= E+ R — (3G — d) (32) 


is used for further discussion of three-dimensional struc- 


tures in this paper. 


Dependent Equations 

In Eq. (32) the question of dependent equations has 
to be answered in a general way before it can be of any 
practical value. The number of dependent equilibrium 
conditions in a set of 6G (or 3G) equations is equal 
to the number of unrestrained degrees of freedom. 
These are the freedoms (independent modes of move- 
ment) of the structure which are not restrained (neither 
rigidly, nor elastically!) by the / elements and K re- 
straints for which columns of coefficients exist in the 
equations. 

It is not uniquely determined which equation, out 
of a group, is dependent. Dependent equilibrium 
equations can be recognized as those corresponding to 
such fictitious supports (expressed in the coordinate 
directions), as would uniquely and sufficiently con- 
strain the degrees of freedom of the model structure. 

Thus, if the structure is sufficiently restrained against 
any rigid-body motion and if it is not a mechanism, no 
dependent equilibrium equations exist. If a structure 
(nonmechanism) is free floating; i.e., considered un- 
supported and unrestrained because it is loaded by a 
group of forces in equilibrium, the equation system has 
six dependent equations corresponding to the six rigid- 
body movements which are not restrained. 

If a structure (or particularly a part structure) is a 
mechanism, it is said to have secondary degrees of 
freedom, one for each independent mode of unrestrained 


motion. An example of such secondary degrees of 
SYMBOLS 
Orem RESTRAINT () \ f . 
O—< RIGID BODY FREEDOM 6, . . 
O—« SECONDARY FREEDOM ant R=6 d=0 
\ J 02 34+6-2°451 
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— (b) G4 
- 7\ 4 
/ | » J R=5 dO 
& yy, A 
4 (¢) t } 
a aT 4 i Fs 
. B7 b4 Pi 
a Q Ps 
6 FICTITIOUS SUPPORTS FOR J R58, del 
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FREEDOMS 


Fic. 3. Degrees of freedom. 
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SYMMETRICAL CASE ANTISYMMETRICAL CASE 


RESTRAINTS AT SYMMETRY PLANE, 
PROVIDED CORRESPONDING LOAD- 
CARRYING ABILITY EXISTS 





ANY PANEL ACROSS THE STRINGERS, PANELS,& BENDING- 
SYMMETRY PLANE IS CONSIDERED ELEMENTS IN THE PLANE OF 
NON EXISTENT SYMMETRY ARE CONSIDERED 

NON EXISTENT 


Fic. 4. Boundary conditions at plane of symmetry 


freedom occurs in a two-dimensional structure. One 
out of every three equations in each gridpoint is de- 
pendent upon the others (and thus is to be eliminated). 
Correspondingly, the equilibrium equations available 
for a two-dimensional structure are Cy) = 2G, and the 
structure has only three rigid-body freedoms. 

A simple case is demonstrated by Fig. 3(a, b,c). It 
is seen from Eq. (32) that dependent equations increase 
the degree of redundancy. Forces in the direction of 
the degrees of freedom (or the mode of mechanism) 
cannot be reacted by the structure. Self-equilibrated 
systems are not impaired by degrees of freedom. 

In a program which automatically formulates three 
(or six in general) equilibrium equations per gridpoint 
of a specified structure, the dependent equations must 
be eliminated from the set to obtain a unique numerical 
solution. Dependent equations must be automatically 
satisfied by the solution; this serves as check for the 
dependency of eliminated equations. Any residual 
from such an eliminated equilibrium equation would 
physically correspond to a support force in the direc- 
tion and at the gridpoint of the eliminated equation. 
This opens up the possibility to express independent 
restraining support forces by the elimination of corre- 
sponding equilibrium equations. In this way, the re- 
straints indicated in Fig. 3 may be expressed by elimina- 
tion of the corresponding equilibrium equations from 


the set. 


Symmetrical Structures 

Symmetrical structures are frequently analyzed by 
using only one half of the structure with the proper 
boundary-conditions on the symmetry plane. Sym- 
metrical and antisymmetrical cases must be handled 
separately. These boundary conditions (Fig. 4) are: 

In the Symmetrical Case—Possible restraints are 
forces perpendicular to the symmetry plane, and mo- 
ment vectors in the symmetry plane. Any panel across 
the symmetry plane must be eliminated because it 
is always unloaded in symmetrical load cases. 

In the Antisymmetrical Case—Two independent force 
restraints in the symmetry plane and a moment re- 
straint vector perpendicular to it may occur. Any 
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Fis. 5. Flat stiffened panel: (a) model of stiffened panel E = 
120, n = 25; (b) determinate structure E = 120 — 25 = 95: 
(c) self-equilibrated internal system used around each internal 
point @ in illustration (a), E = 16. 


elements and bending moment vectors in the symmetry 
plane must be eliminated because they are unloaded 
in the antisymmetrical case. 

All force restraints for symmetry-plane boundary 
conditions can be expressed by eliminated equations 
if the symmetry plane is parallel to one of the coordi- 


nate planes. 


Flexibilities 

The elements are conveniently subdivided into sev- 
eral groups (e.g., wing stringers, rib-caps, fuselage 
stringers, etc.). Frequently, it can be assumed that 
element forces in one element group do not affect the 
element deflections in another. Thus, no cross coupling 
occurs in the matrix multiplications to compute the 
work of deformation [Eq. (14)]. 

The element flexibility D,, of the kth element due 
to the unit load in the pth element is, in general, com- 
puted by the integration JS L,.L,dS, where L;, is the 
unit-load distribution in the kth element (according 
to its definition), and L,dS are the elementary deflec- 
tions in the kth element due to a unit-load distribution 
L, in the pth element. 

The assumption that the deformation of each ele- 
ment is only the result of loads in that element proper 
is, of course, not always true. In each particular 
case, the adequacy of this assumption should be care- 
fully considered. In particular, the incorporation of 
such items as Poisson-ratio influence and swept shear- 
panel effects are important refinements where highly 
oblique panels occur in models of swept or delta wings. 

An oblique panel, loaded by tangential shear forces, 
is under combined elementary shear and normal stress 
and, thus, experiences extension as well as shear de- 
formation. To obtain the correct sweep effect, the 
necessary off-diagonal elements must be provided in 
the flexibility matrix (D), or the corresponding cross- 
coupling matrices D,, [Eq. (14)|. These must express 
the axial deflection of the boundary stringer due to 
tangential shear in a swept panel element, and the 
shear deflection of the panel resulting from axial forces 
in the stringer. The stress-strain relationship in an 
oblique panel has been given by Kappus.’ The oblique 
stresses o;, o,, Tr, and strains e,, ¢, and yz, are given 
in the oblique coordinate system x, y. The sweep 
angle @ is the complement of the angle xy (if xy and the 
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panel are orthogonal: 6 = 0). 
Oblique axial strain: 


| 


" (a 
E-cos 6 


& = 2 — Moy + 2-sin 0-7,,) (27) 


Oblique shear strain: 


2(1 + pv) , Oz + Oy 
VY = E Try + Sin 6 - (28) 


where v = Poisson’s ratio 
w= v— (1+ p) sin’ 6 


A detailed proposal for the handling of both Poisson's 
ratio and oblique panel effects in redundant analysis is 
being prepared for presentation in a future paper. 


(5) Zero Systems and Redundant Internal 
Systems 


In an v times redundant structure, a large number of 
load distributions exist which satisfy equilibrium with 
external loads. Any one of these can be used as a zero 
system [Ly in Eqs. (1) and (4)]. It is not necessary 
to compute the Ly system from a single complete de- 
terminate structure. Rather, it is often convenient to 
choose arbitrarily several determinate load paths, each 
with its own share of the external loads. This results 
in improved accuracy in the solution of the equilibrium 
equations, or permits larger structure to be handled 
within a given program limitation. 

When the external loads are zero, a large number of 
self-equilibrated internal load systems (redundant 
systems) exist. Linearly independent systems, how- 
ever, occur only in sets of exactly systems. Any 
such set of n-independent self-equilibrated systems 
can be used as redundant systems L, in the analysis 
{Eqs. (1) and (4)]. 

One possible set of redundant systems can be com- 
puted from the equilibrium system of the complete 
determinate structure, by reinserting a unit-element 
force alternately into each of the released ‘‘redundant’”’ 
elements. A more general approach is possible which 
has certain engineering and computational advantages. 

Self-equilibrated internal systems can almost always 
be established within a rather small number of elements, 


(a) 





Fic. 6. Internal system for shell structure. (a) eight-panel 
system, for cylinder or conical shells with parallel frames; (b) 
internal system usable on any shape shell structure 
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largely independent of the total size of the model- 
structure. The emphasis in this concept is on the most 
compact arrangement of elements (let this be called a 
substructure) within which a self-equilibrated system 
is possible. The substructure must be one-time re- 
dundant for this purpose. One of its element forces 
is then assigned an arbitrary value (unity, for con- 
venience) for the solution of its member loads. The 
concept of m-redundant elements is replaced by a con- 
cept of n-independent internal systems. 

Eq. (32) for the redundancy m, and the degrees of 
freedom of the substructure are here quite important 
for the formulation and solution of the equation systems 
for the zero systems Ly and internal systems L;. Even 
though the unrestrained degrees of freedom make it im- 
possible to carry corresponding external loads, they do 
not impair the possibility of a self-equilibrated system 
within the substructure. 

Although a larger number of sets of equilibrium condi- 
tions must be formed (one for each internal redundant 
system), certain advantages can be achieved with the 
small systems: 

(a) The number of unknown element loads in each 
equilibrium system is considerably reduced, improving 
computing time and numerical accuracy. 

(b) The smaller size of the internal systems results in 
a larger number of zero terms w;; [Eqs. (4) and (5)]| be- 
cause a number of the systems L; and L; have no ele- 
ments in common. Thus, the coefficient matrix [Eqs. 
(5) and (12)] contains a larger number of zeros which 
improves the accuracy of the solution of continuity 
conditions. 

(c) By decreasing the number of equilibrium equa- 
tions to be set up and solved simultaneously, the over- 
all size of the problem which potentially can be handled 
with given facilities is considerably increased. 

Within the frame of this paper only a few simple 
examples can be given to illustrate the point. Such 
self-equilibrated systems have also been mentioned 
by Argyris.‘ For the analysis of a flat, stiffened panel, 
only one type of internal system is needed (Fig. 5). 

For a cylindrical or conical shell structure, the feasi- 
bility of the internal system shown in Fig. 6(a) (the 8- 
panel system) is easily visualized. If the stringers are 
generatrices of a cylindrical or conical surface and the 
frame segments are parallel, then an open shell is free 
to warp (or twist). Hence, the substructure of this 
system has one secondary degree of_freedom. The 
redundancy of the substructure must be ” = 1 to permit 
a self-equilibrated system: 


Redundancy count: 


22 axial elements, 9 beam elements, and § shear panels. 
Total E = 39. 


Gridpoints G = 15. 
Six rigid-body freedoms plus | secondary freedom: 
d=7 


Thus, n = 39 — (3-15 — 7) = 1 





6 SPANWISE 6 CHORDWISE 
WARPING SYSTEMS 


MODEL, n=2I 


2 AT ROOT 
IN ONE SURFACE ONLY 


3 ROOT WARPING 4 PANEL SYSTEMS 


SYSTEMS 


Fic. 7. Flat multicell wing 


This eight-panel system, which was also used in a 
recent publication® can only be used under these special 
conditions. It cannot be used in noncylindrical or non- 
conical parts, which are not free to warp, because the 
system is then not exactly self-equilibrated. To be 
sure, this twisting rigidity of a general shell surface is 
usually too small for engineering exploitation, but it is 
of importance for the numerical computation. 

Thus, in such areas of general surface (for instance 
in cockpit and tail areas) the substructure with no 
warping freedom must be extended circumferentially 
to a 10-panel system [dashed lines in Fig. 6(a)|. Ona 
cylindrical surface, the 10-panel system should not be 
used because the equilibrium equations will then con- 
tain one dependent equation, leading to a nonunique 
solution. The system, shown in Fig. 6(b) can be used 
in both cases; however, it is larger. The use of one 
fully closed bay which is rigid in twisting in any geo- 
metrical condition, makes the foregoing questions im- 
material. 

Fig. 7 shows the simple set of systems used on a 
multicell wing. The 4-panel system in Fig. 7 is ap- 
plicable in this form only on a flat surface. In the 
presence of (the usual) kick loads along spars and ribs, 
the same system has to be provided with the necessary 
elements from ribs and spars to support these force 
components and make the substructure one degree 
redundant as a three-dimensional space structure. 

As the few examples indicate, this general concept of 
internal systems lends itself to a certain standardization, 
permitting the catalogue listing of standard systems 
for certain types of problems. In complex cases, one 
can frequently recognize the basic problems and, subse- 
quently, modify the standard systems to fit the case 
at hand. 


Cutout Analysis 


Cutouts are a frequent problem in aircraft shell 
structures. Three ways of handling this problem have 
been used successfully. These are briefly discussed: 

(1) Conventional Analysis—The cutout is considered 
as it exists, as an actual discontinuity in the structure. 
The type of internal system, which this requires, must 
encompass the cutout in a frame-like substructure, 
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For plane cutouts, these can be quite simple, as illus- 
trated in the fail-safe crack problem in Fig. 10-V. For 
curved shells, these special systems may become rather 
intricate. 

(2) Infinite Flexibility Method—The structural model 
is made continuous through the cutout, preparing the 
analysis with elements and redundant internal systems 
for the continuous structure. The added elements 
are made infinitely flexible, the load in them becomes 
zero, and the cutout has been effected. However, 
since the actual cutout structure has a lower redundancy 
(by An) than the continuous structure without the 
cutout, An internal systems become dependent, which 
interferes with the numerical solution. Experiment 
has shown that a flexibility increase by a maximum 
factor 10% still leads to a well-behaving continuity 
matrix [equation (12)]. Only An correctly chosen 
elements must be modified (Fig. 13); the remaining ele- 
ments in the cutout then carry almost no load, as well 
as the modified elements. 

(3) Negative Flexibility Method (Fig. 13) 
sary elements can be canceled by providing for each 
a by-pass element of equal negative flexibility. This 
requires An additional, though extremely simple, in- 
Taking, for example, a bar element 
deflection of conventional 


The neces- 


ternal systems. 
and its by-pass (Fig. 13): 
element 

6, = N,(//EA) 
deflection of by-pass element 


5, = No-[—(I/EA)] 


Continuity requires 6; = 62, therefore VN; = —N2. The 
net load transfer through the two parallel elements is 
thus zero, effecting the cut. Again only An elements 
may be canceled in this way for each cutout; a larger 
number of “‘cuts’’ makes the problem physically over- 
determined. 


(6) Program 


A schematic flow diagram of a program mechanized 
for the IBM-704 electronic computer is shown in Fig. 8 
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Fic. 8. Flow diagram of analysis. 
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Fic. 9. Fail-safe test panel. 


It is built around the matrix formulation given in Eqs. 
(10) through (14). For most matrix operations a 
sparse matrix routine was used, which uses only non- 
zero elements. All subroutines for large matrix opera- 
tions and large-order simultaneous equations (depend- 
ing on size of the magnetic core) make use of a con- 
siderable amount of tape operation for storage of inter- 
mediate results. Particularly in the matrix multipli- 
cations to formulate the energy term, the present facili- 
ties of IBM-704 with nine tape units are taxed to the 
utmost. Serving the matrix operations of Eqs. (i0) 
through (14) is a program for the formulation and solu- 
tion of equilibrium conditions for the ” internal systems 
and for the zero-load systems, followed by a sub- 
routine to form the matrices from the result. 

The equilibrium program works from the following 
basic data (Fig. 8). 

(a) The serial specification of the gridpoint coordi- 
nates for the structural model: 


Uf, B, Ch... 14, B, C), 


(b) The serial specification of the elements in m sub- 
groups, by listing the gridpoint numbers which each 
element connects (Fig. 2). 


Axial: (G\G2), 
Bending: (G,G2G3;), 
Shear: (G1G2G3G4), 


(c) The specification of the equilibrium systems by 
serial listing of the element numbers which form the 
substructure of each internal system or zero system. 
The equilibrium equations to be deleted are specified 
by gridpoint number and direction. 

(d) The external applied loads A are specified by 
listing the gridpoint number followed by the three 
load components: 


G,, Ka, Kp, Ke 


Flexibility Specification 
The flexibilities are specified in matrix form from the 
engineering assessment of the properties of the struc- 


tural model. This has been touched upon in previous 
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discussions on the elements of the structural model. 
Proper assessment of the flexibilities of elements in 
the model is of considerable importance in the success 
of any large-scale program. Careful consideration 
must be given to the flexibilities of concentrated rivet 
or bolt connections and of buckled panels. Time and 
space, however, prohibit further detailed pursuit of 
this subject in this paper. 


Checks Provided 

(1) Equilibrium—The equations which were deleted 
from the system of equilibrium equations are auto- 
matically checked with the solution. The residuals are 
printed as well as the solution (element forces) proper. 
From dependent equations, the residual should be zero 
(i.e., small compared with the solution of element 
forces). 

(2) Singularity—If any internal system is linearly 
dependent, the diagonal term of the triangular forward 
reduction of the continuity equations in the corre- 
sponding column will be zero. Zero, however, is a very 
relative matter in numerical computations and, there- 
fore, a comparison with the corresponding diagonal 
term of the unreduced coefficient matrix V is necessary. 
In computer routines which use eight significant figures 
in floating decimal computations, any number which is 
10° lower than another is essentially a zero because it 
consists of mere round-off errors. Depending on the 
size of the matrix V’, lower powers of ten have to be 
treated with suspicion. Thus, the check-routine com- 
putes and provides in the output the ratio: 


diagonal value of V~! 


diagonal value of V 
(3) Iteration—Using the resulting loads (/) of a com- 
pleted analysis as a new zero system in the computation 
of Vy, Eq. (12) provides a direct continuity check: 


Vn = fDF=0 


> 











X* = —V—'-Vy = 0 
With X* several powers of 10 lower than the first set 
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Fic. 10. Typical internal systems 
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Fic. 11. Comparison of strains—panel No. 1 


of solutions, X, a measure of the quality of the coeffi- 
cient matrix and, therefore, of the solution, is obtained. 


(7) Practical Problems 


A few examples of large-scale redundant analysis 
will be briefly discussed. They were selected to repre- 
sent the type of problems which were analyzed in recent 
years with a general program on the IBM-704 elec- 
tronic computer. The program incorporates the fea- 


tures discussed in this paper. 


Fail-Safe Panels 

One scheme for achieving ‘‘fail—safe’’ structure in a 
wing tension surface makes use of a number of longi- 
tudinal strips of integrally stiffened skin planks to iso- 
late and contain a crack. At ultimate strength, a zone 
of plastically yielding material exists in the panel ad- 
jacent to the point of the crack. Attachments near 
the crack are also highly loaded. Static testing to 
account for the many variables is impractical. An 
investigation was made to assess the capabilities of the 
redundant analysis program in handling the nonlineari- 
ties involved in predicting ultimate strength. Six 
static test panels had been run to study this problem 
experimentally. Fig. 9 illustrates a typical test panel 
along with a schematic diagram of one lumped-element 
model analyzed. Various models were tried ranging 
from 59 to 204 redundancies. The five basic types of 
internal systems required are illustrated in Fig. 10. 
An elastic analysis was first completed from which ex- 
cessively loaded stringer elements and rivet “springs” 
were selected. The overloaded elements were handled 
by the simplified method under chapter (2b), with a 
nonreadjusted “tangent line’ approximating the load- 
deflection curve of the critical elements in the plastic 
range. Modified flexibilities were introduced for these 
selected elements and, in some cases, as many as five 
trials were run. These demonstrated conclusively 
that, for the few plastic elements involved, the second 
and certainly the third trial sufficed. 

The computed strains are compared with the meas- 
ured strains for the test panel No. 1 in Fig. 11. Failure 
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occurred in the side panel through the first loaded rivet 
hole next to the crack. The computed rivet -loads 
were comfortably below the ultimate strength of these 
attachments. 

A similar structural model was devised for a repre- 
sentative panel of a box beam test in which rivet failure 
occurred. The degree of redundancy was 50 for this 
analysis; internal systems were chosen among the five 
basic types in Fig. 10. 

The computed rivet load distribution is shown in 
Fig. 12. The first rivet load is indicated critically close 
to the ultimate strength for this type of attachment; 
the side panel, while locally yielded, was not quite 
critical. About four feet of the rivet seam unstitched 
while the adjacent panel material was not visibly 
damaged. 

The redundant analyses of these tests proved straight- 
forward and, with the simple internal systems devised, 
production stress analysis of many variants can be 
accomplished with relative ease. Structural models of 
50 to 60 redundants appear to represent the physical 
problem adequately. Buckling of the skin along the 
crack was accounted for in one panel analysis by in- 
creasing the flexibilities of the affected lateral elements. 
Shear loads combined with tension may be handled, 
although in this case separation into symmetrical and 
antisymmetrical parts is not possible since superposi- 
tion does not hold in the plastic range. 


Fighter Wing-Fuselage Analysis 

(1) The Model—The wing structure is composed of 
sculptured surfaces attached to 15 spars and 6 ribs. 
The wing is attached to five main fuselage rings by 
five highly tapered fittings. The rings are the carry- 
through structure. This structure is represented by 
the model shown in Fig. 14. It includes taper in wing 
thickness and planform, sweepback, and the discon- 
tinuity of the front and rear spars at the side of the 
fuselage. The analysis was performed for symmetrical 
wing loading. 

The complete model contains 491 structural elements, 
138 gridpoints and 33 fuselage supports. From Eq. 
(32) the redundancy is 110. Of these, 69 are in the 
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Fic. 12. Rivet loads—box beam. 
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Fic. 13. Cutouts with negative flexibility. 


wing, 13 in wing-fuselage joint, and 28 in the fuselage. 
The wing internal systems are essentially the types 
shown in Fig. 7. 

(2) Zero Systems—Two determinate systems (Fig. 
15) were used to compute the zero systems. They 
differ only in the part of the wing structure used. 
Each system is essentially a two-spar structure (the 
center spar web and lower surface panels are omitted). 

(3) Test Comparison—The ratio of the local stress 
to the (/ -z//) bending stress is the basis of comparison. 
The chordwise distribution of this ratio for two concen- 
trated unit-load positions (Fig. 16) shows the effect of 
torsion bending. The dips in the calculated curves 
reflect the discontinuity of the forward and aft spars 
at the side of the fuselage (Fig. 14). 

The wing deflections along the 15 and 70% chord are 
shown in Fig. 17 for two loading conditions. As in 
Fig. 16, there is very little torsion effect for the forward 
load, but appreciable twisting along the span for the aft 


load. 


Transport Wing-Fuselage Intersection 

The shell intersection of the wing-fuselage is one of 
the most complex problems of internal load distribu- 
tion problems. The general redundant analysis pro- 
gram was used to analyze the model of a current trans- 
port-type wing-fuselage intersection illustrated sche- 
matically in Figs. 18 and 19. The wing outboard and 
the fuselage portions fore and aft of the center section 
extend slightly more than one wing-chord and one fuse- 
lage diameter, respectively, to insure adequate dissipa- 
tion of the localized effects of the arbitrary load condi- 
tions at the end stations. Taking advantage of the 
vertical symmetry plane, the analysis was made in two 
independent stages for symmetrical and antisymmet- 
rical parts of the load cases. Two unsymmetrical 
and five symmetrical design load conditions, and 
thirty unit-load cases were computed. The symmetri- 
cal and antisymmetrical load separation and the grid- 
point load breakdown into the three components were 
accomplished on a special FORTRAN program (FORmula 
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Fic. 20. Partitioned zero-loads structure antisymmetrical cases. 


TRANslation—an engineer’s program-coding system 
for the IBM-704 computer). The output was a 
punched card deck directly usable in the equilibrium 
equation program for the zero systems. Three separate 
zero systems substructures were devised as shown 
schematically in Fig. 20. This was necessitated by the 
current limit of 240 elements in the program forming 
the equilibrium equations. Portions of the structure 
were omitted where unnecessary for the type of loads 
to be applied; the remaining substructure was, of course, 
determinate. 

In addition to the overall redundancy count [by 
Eq. (32) |, a careful stepwise redundant count was made 
in both the wing and the fuselage. This is most useful 
for the development and location of internal systems. 
The redundancy count summary is as follows: 


Number of Symmetrical Antisymmetrical 
Elements 784 753 
Gridpoints 237 237 
Redundants 178 162 


Many internal systems were of the types shown in Figs. 
6 and 7, modified in various ways according to condi- 
tions. Some special systems are shown in Figs. 21 and 
22. 

Inasmuch as the wing was little affected by the 
boundary conditions of symmetry and antisymmetry, 
approximately a third of the internal systems were 
common to both analyses. Most of the fuselage and 
the wing root required different systems. All grid- 
point, element, and flexibility specifications were 
identical for both analyses. 

Experience with this problem pointed up strongly 
the importance of obtaining the earliest possible auto- 
matic checks on the dependency of the eliminated 
equilibrium equations, and on the independency of the 


internal systems. 


Conclusions 


Some problems in the practical application of re- 
dundant analysis procedures to large-scale complex 
structure have been discussed. These arise mainly 
from the necessity of mathematical exactness in pre- 
senting a uniquely solvable set of equations through 
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the mechanized matrix subroutines which handle the 
numerous computations on automatic high-speed elec- 
tronic equipment. 

The degree of redundancy of the structural model 
must be most carefully analyzed, considering especially 
any secondary degrees of freedom which may exist. 

The dependencies arising in sets of equilibrium equa- 
tions must be carefully eliminated before solution of the 
equations. These later form checking features. 

Emphasis is placed on a concept of providing a set of 
n-independent minimum-sized self-equilibrated internal 
load systems replacing the concept of m-redundant ele- 
ments. Engineering and computational advantages 
arise, particularly in large-scale applications. 

A number of examples illustrate special features 
and demonstrate the versatility of a general approach 
to the large-scale analysis of complex redundant 


structure. 


Appendix (1)—Coefficients of the Equilibrium 
Equations 

Reaction forces in the gridpoints of the three types 

of elements (Fig. 2) are expressed in vector notation. 

Numerical computation follows these formulas in terms 


of the vector components %, y, 2. 


P, = reaction vector in gridpoint G, due to an ele- 
ment force 

g, = locus vector of gridpoint G, 

lx = length of vector g, = g. — 2, le = \En 


Continued on page 614) 





Fic. 21. Wing-fuselage joint systems. Special warping systems 
(modifications of Fig. 7), antisymmetrical case 
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Fic. 22. Special internal systems. 
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On the Thrust Hvpothesis for the Jet Flap 
Including Jet-Mixing Effects’ 


K. T. YEN* 
Rensselaer Polytechnic Institute 


Summary 


This paper is concerned with the thrust generated by a jet 
flap. It is shown that a “‘linear’’ thrust hypothesis can be ob 
tained, provided linearized potential flow is assumed. In fact, 
the linearized problem of a jet-flap system is found to be the 
linear combination of a lift problem and a thrust problem. The 
lift problem gives all the lift generated, but it is of interest to note 
that the thrust problem would yield all the thrust developed by 
the jet flap within the limitation of the linearized theory. 

The mixing of the jet flap with the surrounding fluid is analyzed 
by the momentum-integral method. The analysis substantiates 
Stratford’s suggestion’ for obtaining an increase of thrust by 
causing the jet to mix with the main stream in a region of high 
Finally, some approximate formulas, relating the thrust 
The drag of the airfoil section and 


suction 
and the jet angle, are derived. 
other viscous effects are, however, not considered 


Symbols 

C = chord length of the airfoil section 

f(x) mean camber line distribution of the airfoil section 

g(x) center line of the jet flap 

K(x) curvature of the center line of the jet flap 

p static pressure 

Q mass flow rate 

t(x) thickness distribution of the airfoil section 

‘i thrust 

i direct thrust (Eq. 29) 

", 2 velocity components in x and y directions, respec 

ly 

u’,v’ = perturbation velocity components in x and y direc 
tions, respectively 

5,1 velocity components in x; and y, directions, respec- 
tively 

U, Ua defined by Eqs. (41) and (42) 

W velocity magnitude of the potential flow within the 
jet 

W. = (W, + W.)/2 

“3 = Cartesian coordinates 

Xi V3 curvilinear coordinates attached to the jet flap 

a angle of attack 

€ = width of the jet 

¢ = perturbation velocity potential 

p = density 

0 = angle 

60 = jet angle 
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Subscripts 


0 = the jet exit 


LZ = the upper and lower surfaces of the airfoil, or the 
upper and lower edges of the jet of the mixing 
region, respectively 

= the free stream or the flow far downstream 
a, s = the antisymmetrical and symmetrical problems, 


respectively 


(1) Introduction 


_— of the jet flap have been made by many 
workers in recent years.'~* The lift generation of 
the jet flap has been analyzed by Malavard,' and many 
of his numerical results obtained by electrical analogy 
have been well confirmed.':*7 A similar analysis has also 
been made independently by Helmbold.* The problem 
of thrust has been studied theoretically and experi- 
mentally, mainly with reference to the “thrust hy- 
pothesis,’ which assumes that the forward thrust on a 
jet-flapped airfoil is independent of the jet deflection 
angle.‘~7 (It appears that this ‘forward thrust’’ refers 
to the total jet reaction, which is the thrust developed 
by a jet in a stationary field.) The discrepancy in thrust 
between results predicted by this hypothesis and the 
experimental results was attributed by Stratford to the 
mixing of the jet with the surrounding flow. A one- 
dimensional analysis of the mixing effects has been made 
by Stratford.‘ This analysis is naturally limited in its 
scope. A further study of the thrust problem appears 
to be much needed. 

Some effects of the 
formance of a jet-flap system have been considered by 
Kiichemann.* His analysis is based on an empirical 
relation for the “induced velocity component parallel 
to the airfoil chord’’ due to the action of the jet flap 
[Eq. (44) of reference 8]. In the present paper, a 
formal analysis of the thickness problem is presented. 
However, the main purpose of this analysis is to study 
the thrust hypothesis. 

A complete analytical treatment of the general jet- 
flap problem would bea very difficult task. Malavard’s 
analysis! of the lift problem is based on the assumption 
that the flow is irrotational and the airfoil is thin. It is 
reasonable to consider the flow as nonviscous for an 
On the other hand, it will 


airfoil thickness on the per- 


analysis of the lift problem. 
be seen that viscous or mixing effects play an important 
role in the thrust analysis. However, by adopting the 
same assumption, it is possible to simplify the problem 
considerably, so that the “thrust hypothesis’ can be 


obtained as a basis for further discussion and study. 
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Fic. 1. Jet-flap system. 


In this paper it is first shown that, for an incom- 
pressible potential flow, the linearized jet-flap problem is 
a linear combination of two problems: the lift problem 
and the thrust problem. The lift problem is the same as 
that formulated by Malavard,' while the thrust problem 
is actually the thickness problem in the conventional 
terminology of thin airfoil theory. Since it can be 
shown that the solution of the thickness problem would 
yield all the thrust developed by the jet flap, it seems 
appropriate and meaningful to call it the thrust problem 
for the jet flap. In fact, as the jet angle does not enter 
in the thickness problem, a “linear thrust hypothesis’’ 
follows immediately from this analysis. This appears 
to be an original proof of the thrust hypothesis. 

Next, through a consideration of the momentum, a 
formula for the thrust, in integral form, is given which 
proves to be useful in the present study of the effect of 
interaction or mixing between the jet and the surround- 
ing fluid. The momentum-integral method originally 
developed by von Karman,’ for boundary-layer prob- 
lems, is used to analyze the mixing phenomenon which 
is assumed to be governed by the boundary-layer type 
equations. Although the momentum-integral method 
cannot provide an accurate picture of the mixing proc- 
ess, it is believed that certain overall features of the 
problem obtained through this method are valid. This 
appears to be true for the thrust consideration. 

By using the momentum-integral equation thus ob- 
tained, it is possible to show: (1) how the higher order 
terms modify the “linear’’ thrust hypothesis in the 
potential flow problem, and (2) how the thrust varies 
with jet angle—by considering the mixing between the 
jet and its surrounding stream. Although the potential 
flow theory strongly suggests a more general thrust 
hypothesis, the mixing of the jet with its surrounding 
stream seems to be restrictive and acts to reduce the 
thrust when the jet angle is increased. 


(2) Lift and Thrust Problems of a Jet Flap 


Consider a two-dimensional jet-flap system as shown 
in Fig. 1. Let x and y be the Cartesian coordinates 
chosen so that the free stream is along the positive x 
direction. The free-stream conditions are denoted by 





U., Po, po andy... The airfoil is considered to be thin 
and its chord length is c. The trailing edge of the air- 
foil is chosen to be at x = (c/2) and y = 0. The angle 
of attack is denoted by a. 

In the jet-flap system considered, a high-speed jet 
originating from a source within the airfoil is introduced 
at the trailing edge of the airfoil. At the exit, the jet 
width is €, the jet angle with respect to the x axis 
is #, and the magnitude of the initial jet velocity de- 
noted by W, is assumed to be uniform across the jet. 
The injected mass-flow rate of the jet is Qo = po €o Wo. 

At the exit of the jet, the static pressure fy in the jet 
will be considered equal to that of the surrounding 
fluid. This is in accordance with the following criterion 
as deduced by Courant and Friedrichs:'" ‘‘The total 
thrust is a maximum when the (jet) nozzle is cut off 
at such a place that the pressure at the exit rim just 
agrees with the outside pressure.’’ However, Courant 
and Friedrichs were concerned only with the case where 
the nozzle is stationary. In the usual jet-flap system, 
the static pressure py at the exit may be different from 
that of the free stream. Thus, the actual thrust ex- 
perienced by the airfoil may not be the direct thrust 
poeo WV” cos %. This question will be considered in more 
detail in a subsequent part of this paper. 

It is known that, for an irrotational flow field, there 
exists a velocity potential. As in conventional thin- 
airfoil theory, it is convenient to introduce a perturba- 
tion velocity potential ¢(x, y) for the flow (excluding the 
jet flap) such that the perturbation velocity compo- 


nents are 

u’=u— U, = O¢/x), vw =v = Oe/oy) (1) 
The perturbation velocity potential ¢ satisfies the dif- 
ferential equation 


(0?¢/Ox?) + (0?¢/dy”) = 6 


(2) 


The method of linearization is based on the assumption 
that (uw’/U.) and (v’/U.) are small compared to unity. 
It is necessary, therefore, that (e/c) < 1 and the jet 
angle 4 is smali. 

Let y = fi(x) and y = fo(x) be the equations of the 
prescribed upper and lower surfaces of the airfoil sec- 
tion. The upper and lower surfaces of the jet are as- 
sumed to be described by y = gi(x) and y = go(x). 

The boundary conditions satisfied by ¢ are 


o=u'=v' =0 ato (3) 


v’ U..(dfi/dx), —(c/2) <x < (c/2), y=0, (4) 


II 


I 


2 U.. (dfe/dx), —(c/2) <x < (c/2), y=0 (5) 


v’ = U.(dg,/dx), (c/2)< x, y=0, (6) 


Il 


y’ U..(dg2/dx), (c/2)<x, y=QO (7) 


Since the form of the jet is, a priori, unknown, the 
boundary conditions (6) and (7) will not be used in the 
present form. 

In the present analysis, the linearized potential flow 
problem of the jet flap is separated into a symmetrical 
part and an antisymmetrical part (as in the conven- 
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tional thin-airfoil problem). 
Let v(x, vy) = (1/2) [¢e(x, vy) + o(x, —y)] (8) 
¢a(x, vy) = (1/2) [¢(x, vy) — efx, —y)] (9) 
where ¢, and ¢g, are the symmetrical and antisym- 
metrical parts, respectively, of the perturbation ve- 
locity potential. Let both functions ¢, and ¢, satisfy 
the differential equation (2). 
Let us consider first the symmetrical problem. The 
boundary conditions for ¢, are 


eee oe ae 
Oy/y-0 2L\Ov/, a0, Oy/,—0 


7 ] VU d (f f) = U dt (10) 
o> ae — ae 


— (c/2) <x <c/2 


where ¢(x) = (1/2)(/1 — fe) is the thickness distribution 


of the airfoil section, and 


(S*) =, U de 
oy y=0 7 dx’ 


where e(x) = (1/2)(gi — ge) is the thickness distribution 
of the jet. Since the function ¢(x) is unknown, it is de- 
sirable to transform the boundary condition (11) to a 
more convenient form. It should be noted, however, 
that the variation of the thickness distribution of the 
jet is one of the dominating features of the problem. 
Evidently, this becomes important when a suction re- 


Sx (11) 


wis 


gion is present at the exit of the jet. 

If the flow within the jet flap is also assumed to be 
irrotational, Bernoulli's equation can be applied to the 
upper and lower surfaces of the jet flap. Thus, if the 
density of the fluid is assumed to be uniform every- 


where, 
Me Veen ae oo : 
"3 W,? = a W,? = +> (W,? + 2W,u 2) 

p 2 p 2 p 9 

(12) 
0 l 2 Ls l : 
we = 4 Wet 4 (WW, — Wee) 
p 2 p 2 p 2 

(13) 
where W, = (1/2)(W. + We.) and w, = (1/2)(W —- 


Wy»). Since the pressure must be continuous at the 
interfaces between the jet and the surrounding stream, 
the following two relations are valid 


l l 
E +-,"= + 5 (u* + 11°) ~ 
p 2 p 2 
pi Pele ; 
+ -(U?+ 2Uu,) (14) 
p 2 
] a D2 l - 
2 +a = ; + = (ee + oF) ~ 
p 2 p 2 
2 1 
e+ (U? — 2Uu,) (15) 


where U = (1/2)(u, + we) and wg = (1/2)(u%; — wp). 
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Fiom the above four equations, the following relation 
is obtained, to the first order 


W.? = Woe? + (U*? — U.”) (16) 


where Wx? = Wo? + 2[(po — pa)/p]. Since it is 
reasonable to take the mass-flow rate of the jet as 
peW,, then from the continuity condition for the jet 


eW, = «Wo (17) 
After linearization, Eqs. (16) and (17) yield 


de WoU.d __. ; EWU. [(0*¢; 
a w= 6 = 
dx WW Ox” dx a! Ox . Ox? y=0 


(18) 


and the boundary condition (11) can be written in the 


following form: 


(*) = eWoU. . (S*) 
oy y=0 = Wo? Ox? y= . 


where Wo is a known quantity. Thus the symmetrical 
problem is concerned with the solution of the differential 
equation (2) for ¢,(x, y) subject to the boundary con- 
ditions Eqs. (3), (10), and (19). 

The antisymmetrical problem has been formulated 
and solved by Malavard.' This problem is, of course, 
the “‘lift problem’’ for the jet-flap system. Some of 
Malavard’s results are reproduced here for comparison. 


<x (19) 


Nis 


The boundary conditions for ¢,(x, y) at y = O are 


p if F ; 
(Se) = U ad ” c < ‘ < c (20) 
oy 7 = dx 2 2 
O¢a > 1g , 
( ) av. —, “<x (21) 
oy y=0 dx 2 


where /{(x) is the mean camber line of the airfoil, and 
g(x) that of the jet. Since g(x) is unknown, the bound- 
ary condition (21) is transformed to the following form 
by relating the curvature of the jet with the pressure 
difference across the jet: 


(S¢*) = ct, (es ) . Cc es ra. 
Ox /y=0 i \Oxor/,0 2 


2 2€0 1/2 - 2 po Wo €0 
where €, = (1+ —, ;» {= ar 
ct, boat 


~ 


In Malavard’s analysis, the density of the fluid within 
the jet is assumed to be pp, while that outside of the jet 
ps. 

To evaluate the thrust acting on the wing, only the 
symmetrical part of the problem needs to be considered. 
This can be seen from the following analysis. In the 
antisymmetrical part, the jet flap is reduced to a layer 
of fluid with zero thickness while p!V’*e is kept constant. 
This is necessary, as « — 0, in order to support a fixed 
pressure difference across the jet and to maintain the 
mean camber line of the jet in the same form. Thus as 
€ approaches zero, W? becomes infinitely large (~«~') 
and ¢«W approaches zero (~e'/*). Therefore, in the 
antisymmetrical part of the problem no thrust is de- 
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Fic. 2. Control surfaces and coordinate systems. 


veloped, since the thrust is proportional to the value of 
eW far downstream. Accordingly, a “‘linear’’ thrust 
hypothesis is deduced stating that, within the limitation 
of linearized potential flow theory, the thrust developed 
by a jet flap is independent of the jet angle. In this 
hypothesis the nonlinear and viscous effects are, of 
course, omitted. This is shown clearly by a subsequent 
analysis [see Eqs. (46) and (50) J. 

From the above analysis it follows that a jet flap 
problem, in its linearized form, can be conveniently 
separated into a lift problem and a thrust problem. The 
value of the thrust which is found to be independent of 
the jet angle according to the “‘linear’’ thrust hypothesis, 
may be obtained by solving the thrust problem. How- 
ever, in evaluating the thrust, certain difficulties arise, 
because the solutions of the thrust problem have singu- 
larities at the leading edge of the airfoil section. In 
order to avoid this difficulty, the momentum equation in 
integral form can be used without the restriction of 
linear theory (as in the evaluation of profile drag). This 
is considered in the next section. 


(3) The Thrust Developed by a Jet-Flap 
System: Momentum Consideration 


Consider a control surface for the jet flap system, as 
shown in Fig. 2, which consists of two parts: the first 
part is the circuit 1—-2—3—4, and the second part consists 
of S, the outer surface of the airfoil, and the jet exit 
AB taken perpendicular to the jet flow at the exit. The 
control surface 1—-2—3—4 is assumed to be far away from 
the airfoil so that along 1—2 the velocity is essentially 
U.., and along 3-4 the static pressure is p.. (and the 
density is p..). 

Let Q. be the mass-flow rate across 1-4 and 2-3 
considered negative when the flow is outward. By 
considering the continuity condition of the fluid between 
the two parts of the control surface, it is found that 


0. = po f (u — U.)dy — Qo (23) 


The momentum flux across 1-4 and 2-3 is U.Q.. Let 
7, be the thrust acting on the airfoil by the fluid be- 





tween the two parts of the above control surface. Then, 
from the momentum equation in x direction 


T) = pa f u(u — U,)dy + 


Qi(U. — Wo cos 0) — po «cos % (24) 


The additional thrust 7; can be obtained by consider- 
ing the momentum equation for the control surface 
made of I, inside the airfoil, and AB (Fig. 2). 


Then Tir = Qo Wo cos 0 + po € Cos % (25) 


The total thrust 7 is 


- 
T=7,+ Ty = p. | udu — U,)dy + Qi U 

(26) 
From now on, po will be assumed to be equal to p.. = p. 
The thrust formula (26) is valid for viscous or non- 
viscous flows, and is not limited to linear theory. In 
this form, Eq. (26) may not be convenient for practical 
use, because, as in the “wake survey” problem, the 
integral must be evaluated at a large distance from the 
airfoil (see reference 9). In the present study, however, 
Eq. (26) is found to be very useful. 


Theoretical Value of the Thrust for a Jet Flap With ®, = 0 
and No Mixing 

Consider a symmetrical airfoil section at zero angle of 
attack. If the flow is assumed to be incompressible and 
irrotational, the theoretical value of the thrust for a jet 
with zero jet angle can be easily obtained by using Eq. 
(26). From Bernoulli's equation, the velocity IV in the 
jet far downstream of the injection point satisfies 


W? = Wi? + [(2/p)(po — pa)| (27) 


The thrust is 


T = pW.Woe + of W(W — U..)dy 
— (€/2) 
= p(W,7e), ,e = pW..Woe (28) 
since W,e = Woe in accordance with the continuity 
condition. Let the “‘direct’”’ thrust 7 be defined by 
ti = pW oe (29) 


which is actually the total jet reaction developed by a 
jet in a “stationary” stream. Thus, from Eqs. (27), 
(28), and (29) 


T 27a, when po 2 p (30) 


The above relation shows that the “‘actual”’ thrust 7), in 
general, may not be the same as the so-called ‘“‘direct”’ 
thrust. According to reference 4, for an elliptical airfoil 
section at zero angle of attack, the measured thrust for 
a jet flap at zero jet angle is 94 per cent of the direct 
thrust. Of course, the measured actual thrust would 
include all mixing and viscous effects and a negative 
contribution from the drag of the airfoil. The purpose 
of the above analysis is to show that, when a “‘thrust 
hypothesis” holds, the thrust, in general, cannot be the 
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“direct’’ thrust. Instead, the value 7 = pW..Woe 
given by Eq. (28) should be used (with due considera- 
tion of the drag contribution). In the appendix, it is 
shown that the thrust formula Eq. (28) is consistent 
with the linear theory of Section (2). 

From available experimental data,‘~* it is known 
that the measured actual thrust does depend on the jet 
angle. Stratford’ attributed ‘“‘the discrepancy in 
thrust between the idealized theory and the experi- 
mental results . . . to the mixing with the surrounding 
flow of the two-dimensional jet while still in close prox- 
imity to the airfoil.’” This mixing problem is studied in 
Section (5). In the next section, the momentum inte- 
gral equation is derived for the mixing problem. It is 
of interest to note, however, that some nonlinear effects 
of potential flow on the thrust can also be analyzed by 
using the momentum-integral method. 


(4) Momentum-Integral Equation for the 
Jet-Flap Problem 


As stated before, the static pressure of the jet at the 
exit is assumed to be the same as that of the surround- 
ing fluid. In general, a pressure field will be generated 
by the action of the jet. In addition, lift generation of 
the jet flap depends upon the curvature of the jet. 
Hence, mixing of the jet with the surrounding stream will 
take place under pressure gradients in both streamwise 
and transverse directions. For the treatment of the 
mixing problem, it is convenient to introduce curvilinear 
coordinates x; and y; defined such that y; = 0 is a 
streamline chosen as the centerline of the jet [for the 
precise meaning of the centerline, see Eq. (43)]. Let 
and v,; be the velocity components in the x; and y,; 
directions, respectively. From physical reasoning, it 
is expected that far downstream the directions of x; and 
x will be coincident. 

It is assumed that the mixing of the jet with surround- 
ing fluid will be governed by boundary-layer type equa- 
tions. With the curvature A of the jet of the order of 
magnitude of c~!, where c is the chord length of the 


wing,* the equations are!! 


Ou; Ou 1 Op O71 
U, V; — v = (31) 

Ox; Ov p Ox; Ov; 
Op/oyv; = —pKu;? (32) 
(Ou;/Ox;) + (0v,;/Ov;) = O (33) 


where the curvature A is assumed to be a function of x, 


only. 

Let y; and y. denote the upper and lower edges of the 
mixing region, at which wu; = U,(x;), v; = Vi(x;) and 
u; = Uo(x;), v; = Vo(x;), respectively. 


The momentum-integral equation can be obtained 
by integrating Eq. (31) across the mixing region between 
y. and y;. First, from the continuity equation it can 

* Malavard’s analysis! is based on the assumption that K< 
c—!. In fact, Ke is assumed to be of the order of the thickness or 


camber ratio of the wing. 
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be shown that 


[ Ou, j 
v7,= y 3+) 
4 Ox 
r F dy; d ey 
and VY, = U,— - ujdy (35) 
dx, ax; Jo 
y U dy» d rye j 
2 = 2 _ u,ay 36 
: * dx dx ee — 


At the edges of the mixing region, the flow becomes irro- 
tational—i.e., (Ou,/Ox;) = 0. Hence, from Eq. (31): 


(1/p) (Op/0x;),, + U1 (dUi/dx;) = 0 (37) 
(1/p) (Op/Ox;)y, + U2 (dU2/dx;) = 0 (38) 


From Eq. (32) 


Ae a2 ~ ~ (35) 
_—- dy; = —\y; + 
P Jy, OX; “" Ox 
fay, (pas, 
y - dy; 
B Ov; \p Ox;/ ~ 


_ du, _ dU 
= VV l l -_ yol 2 = [ Vv 
s dx; ; dx, Ce 


By making use of the above relations, integration of Eq. 


(Ku;*)dy; (39) 


(31) across the mixing region yields 


d oy : d "0 ; 
| uu; — U;)dy; + u(u; — Us)dy; 
dx , , F és; J. pe 
dl 1 °y 2 dl 2 0 ’ 
(U; — u;)dy; + | (U2, — u,)dy; — 
dx; Jo dx; J, 


. ” Kyujdy, + K (us WM _ ops =) 10 
- - ’ “y —_ o~Vo ( 
x; In lake NM ox , dx 


The above momentum-integral equation can be used for 
mixing problems in the same manner as its counterpart 
is used in boundary-layer problems.° 
Now, let U,=U+U, (41) 
U2= U—- ty (42) 
where l/(x;) is the symmetrical part and w,(x,;) the 
antisymmettical part of the velocity “‘induced”’ by the 
action of the jet flap. This is analogous to the considera- 
tion given in Section (2) where the potential flow prob- 
lem of the jet flap was separated into a symmetrical 
and an antisymmetrical part. Let the center of the jet 


be defined so that 
ey) wy 
{ ujdy; + } ujdy; = 0 (43) 


Then Eq. (40) can be written as 


d °yi : dl “v1 , 
} uj(u; — U)dy; = (U — u,)dy; + 
ax; J, dx; J, 
; dug 
(41 + ye) (Ung) + (1 — Ve2)Ua - 
P m dx; 2 dx; 


) re) ’ 
vj (Ku;*)dy; (44) 
. Ox, 
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The momentum-integral equation (44) is found to be 
useful in the present analysis of the thrust problem 
concerning nonlinear and viscous effects, since the in- 
tegral in Eq. (26) can be evaluated, in principle, from 
Eq. (44) by integration between x; = 0 and x; = o~. 
Thus, a further analysis of the thrust hypothesis can be 
made. 

It is noted that Eq. (44) is also valid when the mixing 
effects are disregarded. By assuming the flow to be 
incompressible, nonviscous, and irrotational, it can be 
shown (see the appendix) from Eqs. (26) and (44) that, 
after omitting some higher order terms, 


d : 
| o + yo) (Uu,g) — 
dx; 


| yj (Ku; )dy; | dx; (45) 
v2 Ox; , 


This checks with Eq. (28) for the case AK = u, = 0. 
For Kc < 1, Malavard! has shown that 


T= pW.Woe + p | 


/70 


Ug ~ Ke(W.?/U.) 


Hence, it follows that, according to potential-flow 
theory, from Eq. (45) 


L = pW..Woe[1 oo O(K je) | (46) 


where A; is a typical curvature of the jet flap and 0(K je) 
denotes terms of order K je or higher. Thus, the linear 


“thrust hypothesis” is valid to the order of A je. 


(5) Effect of the Jet Mixing on the Thrust 


The analysis given above strongly suggests that the 
deviation from the thrust hypothesis is mainly due to 
viscous effects which invalidate the assumption of 
potential flow. It is proposed, therefore, to analyze the 
mixing of the jet with the surrounding stream and its 
effects on the thrust. This problem is, however, a com- 
plicated one, and the present study will make use of 
Eqs. (26) and (44). Obviously, without attempting to 
obtain the solution for u;, U, 1, v2 etc., a complete 
understanding of the phenomenon is not possible. In 
addition certain plausible assumptions have to be in- 
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Some qualitative discussion of the mixing problem 
will be given first. Since for the usual jet-flap system 


[ (u; — U) dy;> 0 


the first term on the right hand side of Eq. (44) is 
either positive or negative in accordance with the con- 
dition dU/dx; S 0. Thus, if the pressure at the exit is 
lower than that of the free stream, dU’/dx; may be ex- 
pected to be negative for all x;._ Then the first term in 
Eq. (44) is positive and is favorable for thrust genera- 
tion. Further, for dU/dx;< Oand u,> O (for all x;), it 
is expected that y; + y2> 0. Thus, the second term is 
also positive, although it is of the order Ae* for Ac < 1. 
This supports the suggestion of Stratford! that a gain 
of the thrust is achieved by allowing the mixing to 
occur where the pressure is less than the free-stream 
pressure. On the other hand, a loss would result for 
mixing which occurs where the pressure is higher. In 
Eq. (44), the third term on the right hand side is, how- 
ever, generally negative. But this term is of the order 
K*e*, while the last term is of the order A e*, where A is 
the curvature and ¢ is the width of the jet flap. The 
contribution from the last two terms is small when 
Ke < 1. (For large x, it is expected that A ~ x™, 
at most, e ~ x"/?, and Ke? ~ 0Oasx— o.) 

Since the derivation of quantitative information from 
Eq. (44) is a formidable task, the following analysis 
will be limited to small 6 and small jet curvature. 

For a jet flap with small jet angle % and small curva- 
ture A, the last three terms in Eq. (44) can be neglected. 
Thus, to find the thrust, only the integral of the first 
term on the right-hand side of Eq. (44) needs to be con- 


sidered, hence 


yy) ; ey dU vy ; 
| u,j(u; — U)dy; =| | (U — u;)dy; | dx; + 
J y2 ij 0 dx; Jy z 
Woeo( Wo — Uo) 


In order to gain some information about the de- 
pendence of the thrust upon the jet angle (without de- 
tailed knowledge of the functions u;, LU’, ete.), it is pro- 
posed to adopt a simplifying assumption, namely, that 
the thrust dependence upon the jet angle is mainly due 
to the action caused by the deflection of the jet flap, 
while UL’ is relatively insensitive to changes in %. This 
appears to be physically reasonable for small 4) and is 
not expected to introduce serious errors for the thrust 
evaluation. (G.I. Taylor has adopted a similar assump- 
tion for a vertical-jet problem. See reference 1.) 


Let vy = g(x) be the centerline of the jet flap, and 
tan 6 = g’(x) (Fig. 3). If the jet angle is small, the 


following relations are valid: 


du dU 


dx; dx 


cos 6 


Ni 


dy; ~ dy cos 8 


Thus the integral on the right-hand side of the preced- 
ing equation can be approximated by 
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THRUST HYPOTHESIS 


so «6aU 1 
(cos? 6) | | (U- ait) dx 
J (c/2 dx J ¥ ‘ | 


and the thrust is given by 


] T, - Q(U. a Uo) + 


F dl € *y1 , — 
p (cos* 6) | (U —u)dy|dx (47) 
é 2 dx Jy 


where 74 = plVo"e is the “direct thrust’”’ and Qp is the 
total mass-flow 1ate of the jet at the exit. Eq. (47) 
shows that the thrust dependence on the jet angle is 
accounted for by the factor cos? 6. Thus, for a sym- 


metrical airfoil section at a = 0 


Pa 6 = T, + Qo( U om Uo) + 
coal | "(U ty|dx (48 
ons ) 1| dx 4! 
of 2) dx fe tel iy ; 
If no mixing is considered, 74)» = pW..Woeo (see the 


appendix). Hence, it is evident that the “‘linear’’ thrust 
hypothesis is implied in the case 6 ~ 0. It should be 
noted, however, that for an arbitrary airfoil section, the 
value of 6 is not necessarily zero when a = 6 = 0. 

A simple approximate formula for 7 can be obtained 
from Eq. (47) by replacing cos? 6 by the value 


[1 + a cos? (6 — ba)|/(1 + a) 


where a and 0 are “shape parameters’ presumably to be 
determined from experimental results. Thus, 


l 
T= 7T74+Q(U. — Us) + xX 
os l+a 
: . dU ny; 
[1 + a cos? (@ — ba)| | | (U — u)dy \dx 
c/2 dx 42 
(49) 
which can be written as 
1+ 20)T, + 
= ( a Cos- =() 
l+a ” 
: (1 7OA)(T, + Q(U Uy)] (50) 
—- oo ’ ( i a ( 3) 
l+a 
where 9 = 0 — ba, and 79-9) = 74,-) as given by Eq. 
(48). 
For thin airfoil sections at small angles of attack the 
contribution from the term Qo(U.. — WU») may be as- 


sumed to be small. If this term is neglected, the follow- 
ing simple relation is obtained: 


l aes 
T= (i+a cos? 8)T» Oi ia 
l+a 
a " oan - 
ae (1 — cos? @)7, (51) 


In a comparison of the above formulas with the 
measured actual thrust data,® it was found that the 
values of thrust obtained from Eq. (51), for example, 
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are too high.* Of course, it is to be kept in mind that 
the measured actual thrust values include the drag con- 
tribution and other viscous effects. The above analysis, 
however, accounts for the effects of mixing between 
the jet and the surrounding stream. 


By comparing the “‘linear’’ thrust hypothesis, equa- 
tion (46) and the results obtained in this section, some 
tentative conclusions concerning the thrust hypothesis 
may be made. It appears that, strictly speaking, the 
only valid hypothesis which can be proved rigorously is 
the ‘‘linear’’ thrust hypothesis. However, if the flow is 
assumed to be nonviscous and irrotational, and the 
curvature of the jet flap is small, the thrust hypothesis 
will be valid for fairly large values of the jet angle. 
However, by considering the mixing between the jet 
and the surrounding stream, it appears that the thrust 
is sensitive to a change in the jet angle. Thus, it is 
evident that the mixing phenomenon plays a very im- 
portant role in the thrust developed by a jet-flap 
system. 


(6) Concluding Remarks 


The results obtained from this study may be sum- 
marized as follows: 

(1) The linearized jet-flap problem can be separated 
into a symmetrical and an antisymmetrical problem. 
The symmetrical problem accounts for all the thrust of 
the jet flap and a “linear” thrust hypothesis 7 = 
pW..Woe follows directly from this analysis. 


(2) A simple integral formula for the thrust de- 
veloped by a jet-flap system is given. In general, the 
“direct thrust’ developed by a jet situated in a station- 
ary field is not necessarily the same as that actually ex- 
perienced by a jet-flap system in a uniform free stream. 
This discussion is based on a simple analysis without 
considering viscous effects. 


(3) In accordance with a higher order analysis, it 
was found that, when potential flow is assumed, the 
thrust may be approximated by 


T = pW.Woel[l + 0(K,«)] 


where A, is a typical curvature of the jet flap. There- 
fore, the thrust hypothesis is limited to linearized jet- 


flap problems. 


(4) Stratford’s suggestion that a gain in thrust can 
be obtained by allowing mixing to occur in a region with 
suction is substantiated by a more general analysis 
carried out in this work. 


(5) By considering the mixing effects of the jet with 
the surrounding fluid, formulas are obtained for the 
thrust dependence on the jet angle (for small jet angles 
and small jet curvatures). 


* However, by fitting the experimental results (Fig. 1 of ref- 
erence 4) with Eq. (50), fairly good agreement can be obtained up 
to 9% =~ 50°. 
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First, a proof will be given to show that the thrust 
formula 7 = pW... Woe is consistent with linearized po- 
tential theory for a symmetrical airfoil section with a 
jet flap at zero jet angle. The addition of Eqs. (12) and 
(13) yields to the first order 

W? + [(pi + pe)/p| = const. 
Similarly, from Eqs. (15) and (16), 
U? + [(pi + pe)/p| = const. 
Subtraction yields 
W? = U* + const. 


Thus, from Eq. (44), by making use of the above rela- 
tion and the condition of no mixing, which can be ex- 


a 
of u,dy; = Qo = const. 
y 


it is found that 7 = pW..Woe. 
Similarly, the proof of Eq. (45), wherein only the w,” 
terms are omitted, can be carried out by using the same 


pressed by 


equations. 


AEROSPACE 


1960 


SCIENCES ATG UST, 


References 


1 Malavard, L., Contribution a l'étude theoretique du soufflage au 
bord de fuite d'un profil d’aile, Note intérieure ONERA 4/1727A, 
December, 1954. (Also AD75662. ) 

? Helmbold, H. B., On the Lift of a Blowing Wing in a Parall. 
Stream, Journal of the Aeronautical Sciences, Vol. 22, No. 5, pp 
341-342, May, 1955. 

3 Spence, D. A., On the Lift of a Blowing Wing in a Parall 
Stream, Journal of the Aeronautical Sciences, Vol. 23, No. 1, pp 
92-94, January, 1956. 

‘Stratford, B. S., Mixing and the Jet Flap, The Aeronautical 
Quarterly, Vol. VII, Part 2, pp. 85-105, May, 1956 

> Dimmock, N. A., An Experimental Introduction to the Jet 
Flap, A.R.C. Current Paper 344 (British), February, 1956 

6 Dimmock, N. A., Some Further Jet Flap Experiments, A.R.C 
Current Paper 345 (British), February, 1956 

7 Williams, J., British Research on the Jet-Flap Scheme, 
Vol. 6, No. 6, pp. 170-176, June, 1958. 
1 Method for Calculating the Pressure Dis 
Royal 


ZEW 


’ Kiichemann, D., . 
tribution over Jet-Flapped Wings, Rept. No. Aero. 2573, 
Aircraft Establishment, May, 1956. (AD 101815.) 

® Schlichting, H., Boundary Layer Theory, McGraw-Hill Co., 
1955 

1 Courant, R., and Friedrichs, k. O., Flow and 
Shock Waves, pp. 392-394, Interscience Publishers, New York, 
1948. 

11 Murphy, J. S., Some Effects of Surface Curvature on Lamina) 
Boundary-Layer Flow, Journal of the Aeronautical Sciences, Vol 
20, No. 5, pp. 338-344, May, 19538 


Sut ersonu 


Analysis of Redundant Structures . . 


(Continued from page 606) 


(1) Bar element G\G: 


4 ee 
ee 


(2) Beam element G,G.G;. Unit bending moment 
vector /: is normal to the plane of G,G.G;; the plane 


must be specified if the three points are in line. 


Si = 8 — B11, Bes = Ss — 
ho = Zi X Bex (vector product) 
h= ks Unit moment vector / is iy divided by its 
absolute value |i 
P, = -— (go X h), Ps = — ra (2s X h) 
P, = —(P, + P) 
(3) Panel element G\GoG;G, 
oe sa foe - es - & 
Diagonal forces: 
Ly: = iin le =~ “gees 
2 2 


Gridpoint reactions: 


3 l 
P, = —P, = t= / _ _ a q(&s a £1) 
13 - 
= ss 1 
P, = —P, vo Las : i 518s —€ 
24 “ 
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Minimum Wing Wave Drag With Volume 
Constraint 


T. STRAND* 


Convair, A Division of General Dynamics Corporation 


Summary 


\ numerical method is developed for calculating the minimum 
a given wing planform and volume using 
The corresponding optimum 


thickness drag for 
linearized supersonic flow theory 


volume distribution is also determined. The results show that 


considerable drag reduction is possible by improved volume 


distribution 


Introduction 


i ies DIRECT calculation of the zero lift pressure drag 
(or wave drag) of given wings using linearized 
supersonic flow theory has occupied many authors. 
A listing of these authors as well as summary charts 
these numerous 
It may be noted 


of the wave drag compiled from 
sources can be found in reference 1. 
that the wings considered, although varying in plan- 
form and section characteristics, all have spanwise 
constant thickness ratio. 
Except for special cases—e.g., wings with wedge sec- 
tions or taper ratio of unity—the wave drag calculation 
is lengthy and tedious and therefore best suited for 
digital computer programming. Some time ago Con- 
vair—San Diego completed an IBM-704 routine for 
the calculation of wave drag of wings with straight 
leading and trailing edges, parabolic-arc airfoil section 
and incorporating the possibility of varying the thick- 
ness ratio linearly along the span. The approach used 
was essentially an adaptation of the numerical method 
given in reference 2 for the calculation of drag due to 
lift. This method was again an extension of a pro- 
posal made by Dr. Ta Li in reference 3. In reference 
| this IBM-704 program was employed to obtain the 
wave drag of a large number of wing planforms with 
parabolic sections at Mach number 3. It was found 
that it was not necessary to perform a systematic in- 
vestigation of the effect of different types of linear thick- 
ness ratios, as the data points correlated perfectly on 
the basis of the effective thickness ratio defined as fol- 


lows: 
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2 


9 20/2 
— } (t/C)? Cdy 
SJ0 . 


(7 C) ase” = (1) 


Here S is the wing area, ) the span, C the local chord, 
and ¢ the maximum thickness. 

This conclusion was further strengthened through 
additional runs by the present author, who also showed 
that if the drag coefficient was a function of the effective 
thickness ratio (a hypothesis which, by the way, could 
not be verified analytically), the minimum wave drag 
for a given volume would be obtained by designing 
wings with the spanwise thickness ratio distribution 
proportional to the local chord. This is then a cri- 
terion based on the results of numerical integration 
only, and is valid for wings having circular are sections 
and spanwise linear taper in thickness ratio. 

The problem of finding the minimum wave drag for a 
given volume without restrictions on the section proper- 
ties and the kind of spanwise thickness taper was still 
left unanswered. The solution to this problem is the 
subject of the present study. 
ment an account will be given of the numerical method 
by which the absolute minimum wave drag of certain 
types of wing planforms with specified total volume was 


In the following develop- 


obtained. 


Mach Box Method 
Drag 


The basic ‘“‘Mach box grid’’ method of computing 
drag has been adequately treated.” However, for 
completeness in presentation and because there are 
some minor changes in approach, the method will be 
briefly reviewed here. 

Let us start by assuming that a given wing planform 
has been selected and that we have an analytic ex- 
pression for the surface ordinate (2) normal to the 
plane of the wing. The wing and the disturbed flow 
region adjacent to it in the plane of the wing is now 
overlaid by a rectangular Mach box grid as shown in 
Fig. 1. The term Mach box originates from the fact 
that all the diagonals of the boxes are Mach lines. 
The downwash within each box is considered to be 
constant and equal to the value in the center of the 
box. The velocity potential at an arbitrary point 


(x,y) can then be written as follows: 
o(x,y) = 


U Oz dé dn 
226) ioe 
Ww j Ox j 19 Vix — &)- B?(y —n)° 


(2) 
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where L’ is the free stream velocity, 02/0x is the down- 
wash angle in the center of the box, and 6 = (/? — 
1)'/?. The area of integration A; is the area of each 
box lying within the forward Mach cone from the 
point (x,y). Next we introduce nondimensional co- 
ordinates by dividing by 5;, the streamwise length of 
each box, and we eliminate 8 from the kernel in the 
following manner: 


xy = x/b,, yn = By/bhy (3) 
& = &/h, m = Bn/b 
Then 
bU Oz 
(X1, V1) = 28 > (=) x 


[—(2 m)|d& dn 
{f, cose 
1j V (x — §)? — (wn — m)? 


Letting vy = x, — and yu = y, — m, it is seen that the 
integral part of the right-hand side of the expression 
above becomes 


2 dv du 
- T Aj V/ 2 a we ( 


The center of the box in question is indicated by the 
coordinates (subscripts) 7, @ (Fig. 1). We note that we 
have three types of boxes. Hence for each type we 
obtain a separate expression from the integration 
namely, the following influence numbers: 


vt 


of Eq. (5); 
Roo =] 


sin! u 1 — V1 — p’\ |u= a 2 
— In b+ (1/2) 
Me i u=l 


is valid for 7 > 0 


2 l 
Ria = (« — 5) xX (6) 
Tv o 
= lu | aS p> Fie 
— In é @ 
u u p> Eas 
2 l 
~2 (p+ 2) 
T “ 
a+(1/2) 


sin"*# _ l~ Wi =»? P= (1/9) 
" ws s n= Bt(1/2) 


y— (1 


The last expression is valid fori <7> 0. Substitution 
of Eq. (5) into (4) yields 
286(%1, V1) Oz o 
Ao = Box, ¥1 a ¥ (R ) (7) 
byl Ox) 5,5 


where the summation is to be performed over all the 
boxes within the disturbed flow region in the forward 
Mach cone from the point (x, 1). 

The drag coefficient is given by 


oS ee (S) 
U sax ox ig 


Cas — 
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where 2 represents the surface ordinate normal to the 


plane of the wing. 
Integration by parts yields 


: 4 ‘f Oz 
i == i a 
Cr U Wre ? Ox ” 
Oz 


ra O72 ) 
J, E % ox dy — ff. g x? dx dy¢ 


(9) 
_ 4 fal Oz =] 
—  UlB p> "=. p» “ee 
b,? 0° | 
B X ? ox?f 


where L.E. and 7.E. indicate the leading and trailing 
edges, respectively. Substitution of Eq. (7) into (9) 
gives finally the following expression for the drag 


coefticient: 


Oz OZ 0227] 
. as 4 - A 10) 
is Ox p» ? ox 1 2X ? 4 ( 


Boundary Condition 


The boundary conditions for the computation of the 


velocity potential [Eq. (7)| are very simple. Since 2 
is the local wing surface ordinate we have 
Oz O02 : 
= on the wing (11) 
Ox z=0 Ox 
Oz ‘ : 9 
= () off the wing (12) 
Ox|,=0 


Thus since the downwash angles are known every- 
where, the problem of finding the drag coefficient of a 
given wing surface is reduced to a fairly straight- 
forward numerical calculation procedure involving the 
determination of the velocity potential at a number of 
spanwise and chordwise discrete points, multiplying 
the potential by the first or second streamwise deriva- 
tive of the local wing surface ordinate and summing 
the products using Eq. (10). 

This is then a brief summary of the numerical method 
reported in reference 2 which provided the basic tool 
on which the subsequent original development is 


based. 


Volume 
Because it is preferable to work with nondimen- 
sional coefficients, we shall define a volume coefficient 


as follows 


C, = V/SC (13) 


where V is the total wing volume, and C the mean 
aerodynamic chord. Notice that since V/S is the 
average wing thickness, C, is a kind of average thick- 
ness ratio. 

The expression for the wing volume coefficient in 
terms of the given wing coordinates is now 
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MINIMUM WING 


2 
C= “ff 2 dx dy (14) 
SG S ‘ 


which in the Mach box summation system becomes 


«SE ey (15) 
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Similarity 

In order to be able to plan ahead for the drag optimiz- 
ation procedure it is advantageous to know the group- 
ing of the basic parameters involved. Hence a simi- 
larity analysis is indicated. It can easily be shown 
that the ratio of the wing surface coordinate to the 
volume coefficient may be written 


z/C, = f(x,y) (16) 


where f(x,y) is a known expression. [For a wing with 
a parabolic are section and thickness ratio proportional 
to the chord, f = —4(x — xz.%.)(x — xr.2.)/C,;. Note 
that in this case (t/C)?.,, = (3/2)(t/C2)EC,. Here ¢, 
C and C, are the local maximum thickness, local chord 
and the wing root chord, respectively. | 

A similarity analysis for planar flow similar to the 
one outlined in reference 6, shows that for the same 


shape function f 
Cp = Cp (C,?/8, B &) (17) 


Now, since the pressure coefficient is of necessity 
directly proportional to C,/8, the drag coefficient must 
be proportional to C,’/8 and we are justified in writing 


BC)/C,? = fn (8 &) (18) 


This is therefore the manner in which the results will 
be presented. It may be noted that a two-dimensional 
wing with circular arc or symmetric double-wedge 
airfoil section has 8Cp/C,* equal to 12 or 16, respec- 


tively. 


Optimization 


Series Expansion 

The wing surface ordinates are next expanded into a 
series, each term of which has the property that the 
ordinate vanishes along the leading and trailing edges. 
A ten-term series, chosen more or less arbitrarily, which 
has this property, is shown below: 


q 


/ —- f 
= >. Q(z) 


oa 


7=0 
where by definition 
% =A 3, = Ax’ |? 
Zz,’ = A y’ . Ze’ >= Ax y’ ’ 
Se = A iy’ i* 3,’ = Ax? 
33" = A y’ ‘ Z,' = Ax? y’ 2 
Ss,’ = Ax Zo’ = Ax? 
Here 2’ = 2/C,, y’ = 2y/b and A = —(x — xz.2.) X 


(x — xr.z.)/C,?.. The a,;are unknown coefficients which 
will be determined by the requirement that the drag 
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for a given wing volume be a minimum. Notice that 
the first term in the expansion represents a wing with 
circular are sections having a spanwise thickness ratio 
distribution proportional to the local chord. Thus we 
would expect from the remarks made in the introduc- 
tion that most of the possible drag reduction would be 
obtained by this term alone. This is indeed the case, 
as will be shown later. 

Each 2,’ above represents by itself a certain wing sur- 
face shape with its own drag coefficient which we may 
call Cp,,. Associated with each term 2,’ there will 
also be a pressure distribution which will interact with 
the angle of attack distribution of every other term in 
the series. Thus the drag coefficient may be written 


y 
Co= > aan, (19) 

i=0 ; 

7=0 


This is thus the total pressure drag coefficient due to 
wing thickness. 

Similarly the total wing volume coefficient is ob- 
tained by summing the volume coefficients of each 


separate wing 2,’, i.e., 


C,= >) ac,, (20) 
+=0 
Here Cp,, and C,, are obtained from Eqs. (10) and (15), 


respectively. 


Minimum Drag With Volume Constraint 

The problem at hand may now be formulated as 
follows: It is required to find the wing surface, ex- 
pressible in terms of the given series expansion, having 
minimum drag coefficient for a given volume coefficient. 
This is a problem from the ordinary theory of maxima 
and minima. The usual approach is now to define a 
new function F which has been chosen with due regard 
to the previous similarity consideration : 


F = BCp + uC, (21) 


Here yu is a constant Lagrangian multiplier. The ex- 
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tremum of Cp is obtained when the function F fulfills 
the system of linear algebraic equations 


OF 
= i i (22) 
Oa ; 
OF 
- = () (23) 
Ou 


By multiplying Eq. (22) by a,, i.e., by multiplying 
the first equation of (22) by dp, the second by a; and 
so on and then summing the resulting set of equations,’ 
it is easy to show that 


BCp = uC,/2 (24) 


For computational purposes it is expedient to choose 
C, equal to unity and then get a solution to Eqs. (22) 
and (23), i.e., to determine the a; and the Lagrangian 
multiplier explicitly. 

It can also be shown from Eq. (22) that yu is directly 
proportional to the volume coefficient C,. Hence we 
may write 


B Cp/C,? = —(1/2)ucp=1 (25) 


The set of a; determined by setting C, equal to unity 
must now be used with the surface ordinates per unit 
C, or 


2’/C, = > Qin 1°21 (26) 
7=0 

It should be mentioned here that in the optimization 
scheme outlined above there is no built-in guarantee 
that all surface ordinates of the wing with minimum 
drag will be positive. The possibility of getting nega- 
tive ordinates was considered at the start of the present 
investigation, but was rejected as being most unlikely 
Negative volumes over part of the plan- 
These have 


to happen. 
form did, however, occur in some cases. 
been discarded. It was learned much later that the 
possibility of getting negative ordinates in optimum 
volume cases had previously been pointed out.’ 


Discussion 


The Mach box grid size required for adequate accu- 
racy was determined first. Figure 2 shows that the 
constant level is for practical purposes reached with 
about 16 boxes in the streamwise direction. This grid 
size was then used in the subsequent runs on the digital 
computer. The results for three different series of 
wing planforms with taper ratio zero are presented in 
Fig. 3. In addition to the minimum wave drag for a 
given volume, also shown are the drag values for the 
same wings with parabolic are sections and spanwise- 
constant maximum thickness ratio, and with maximum 
thickness ratio proportional to the local chord. It is 
seen that considerable drag reduction is theoretically 
possible with the use of improved volume distribution. 
It is also significant that a major part of the possible 
drag reduction is obtainable by use of parabolic arc 
sections and a maximum thickness ratio distribution 
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Fic. 5. Optimum thickness distribution, BR = 5 


proportional to the local chord variations. Since the 
parabolic are is the section for minimum wave drag 
for a given volume in two-dimensional flow, all curves 
tend toward BCp/C,” = 12 in the limit of large BR. 
Sample thickness distributions for minimum wave 
drag are shown in Figs. 4, 5 and 6. Figure 4 indicates 
a feature which is common to all three types of wings 
as 6 A is reduced, namely, an unexpected decrease in 
thickness over a certain part of the wing. For the wing 
shown in Fig. 4 this necking-down of thickness occurs 
around the midspan. As 8 R approaches that corre- 
sponding to sonic leading edges this feature becomes 
more and more pronounced and finally a point is reached 
after which we actually have negative volume over part 
of the wing. The above tendency is not yet apparent 
in Figs. 5 and 6, as 8 MR is not low enough. It was 
found, however, that all minimum drag wings with 
sonic leading edges exhibited negative volumes over 
part of the wing area. This is of course a physical im- 
possibility and is the reason why the curves in Fig. 3 
stop short of the 8 M corresponding to wings with sonic 
leading edges. The reason for the occurrence of nega- 
tive ordinates is not clear, but it must obviously result 
in a better overall area variation of the equivalent 
bodies of revolution. In any case, it can be concluded 
that wave drag minimization by the ten-term series 
method used here is feasible for wings with supersonic 
leading edges only above a certain lower limit in 6 &. 
A calculation of the optimum wing ordinates using the 
two first terms only in the series expansion also showed 
part negative volumes for wings with sonic leading 


edges. 


Conclusion 


A numerical method has been developed for calcu- 
lating the minimum thickness drag of supersonic wings 


having a given total volume. The method is applicable 
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Fic. 6. Optimum thickness distribution, 8 R = 7 


to wings of supersonic leading edges above a certain 
8 MR. The results show that considerable drag reduc- 
tion is possible by improved volume distribution. The 
major part of the drag decrease can be obtained by using 
parabolic are sections with the maximum spanwise 
thickness ratio proportional to the local chord variation. 
Hence, if volume is a criterion the root thickness ratio 
of these wings should be much higher relative to the tip 
thickness ratio than is common practice at subsonic 
speeds. 

For wings with sonic and subsonic leading edges the 
optimization procedure failed in that negative wing 
ordinates were indicated over part of the wing. 
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The Thickness of a Melting Ablation-Type 
Heat Shield 


Ernst Wilhelm Adams 

Aeroballistics Laboratory, Army Ballistic Missile Agency 
Redstone Arsenal, Ala 

March 28, 1960 


—_—— FOLLOWING DISCUSSION shall deal with some material 
aspects of the ablation-type heat protection for a certain 
assumed IRBM with a homogeneous, opaque glass shield which 
is allowed to melt and evaporate. For the vicinity of the 
stagnation point of this shield, a new numerical procedure has 
been derived recently (see references 1 and 2) which yields 
highly accurate transient solutions of the pertinent system of 
partial differential equations. The solutions consist of the 
temperature 7(t,z)[°K.] and, for the liquid film, the axial 
velocity component v(t,z) [mm./sec.] as functions of time 
t [sec.] and axial distance z[{mm.] normal to the surface. The 
reference system is indicated in Fig. 1. 

The calculation procedure includes the assumption of a semi- 
infinite shield, 0 < zs < o, with the boundary condition 
lim 7(t, 3) T) = 300[°K.]. If a thermally insulated inner 
edge of the shield is placed at a station where the ‘‘semi-infinite”’ 
solution indicates, at a certain time instant, say, 380[°K.], then 
the true temperature at this station is about 460[°K. ]. 

The necessary weight of the shield is governed by (1) the time 


ly 
integral s = -f 
0 


rate vo(t)[{mm./sec.] from the beginning of the re-entry at 
t = O[sec.] to the impact of the vehicle at ¢ = ty[sec.]; and (2) 


Vo(t) dt [mm.] extended over the melting 


the maximum thermal penetration across the shield, which 
evidently occurs at impact time, and is measured as the distance 
from the surface to the station 6 = b(t;)[mm.] at which the 





Contributors Please Note 


Readers’ Forum items must be confined to the 
equivalent of one page in the Journal. Contribu- 
tions that exceed this limit will be returned to the 
authors for condensation and rewrite. To avoid incon- 
venience and delay, please adhere to the following 
specifications: 

(1) Five, double-spaced, typewritten, manu- 
script pages (82 by 11 in.), including 
wide margins, formulas, and headings, 
equal one printed page. 

(2) For every illustration, deduct at least 
one-half or as much as one manuscript 
page, depending on the size of the 
illustration. 











a 








n 
h 


Thee me 


a 


ee 





solution reaches ZT = 880[°K.]. The minimum 
thickness of the shield is given by the sum (s + }){mm. ]. 

For the IRBM re-entry under discussion, solutions have been 
calculated on an IBM 704 computer for a wide variety of the 
relevant material properties including the measured viscosity 
functions 


necessary 


functions w(7) [kg. sec./m.?] and vapor pressure 


p(T) |kg./m.*] of Pyrex glass and Fused Silica. For Pyrex 
glass, the level of both w(7) and p,(7T) has been varied by use 
of nondimensional factors u*{— | and A[— ], 


w*Xupyrex(7T) and p,(T) AX Po. pyrex(T 


For the sake of simplicity, non-temperature-dependent but varied 
constants have been used for the other relevant material proper- 
ties, i.e., the thermal conductivity k [keal./m. °K. sec.], the 
specific heat c, [keal./kg. °K.], the specific weight + {kg./m.*], 


and the emissivity constant e{[—] of the surface of the opaque 


u(T) 


shield. 

Fig. 2 presents some results of this material investigation 
Each calculated point, indicated by a circle, triangle, or square, 
stands for the data b |mm.| or (s + b)[mm.] of one and, in some 
cases, more than one of the 35 investigated glass shields with 
different material properties. These points allow to draw one 
(solid) line for (s + 6)[mm.] and several (dashed) lines, with 
u* as curve parameter, for b{mm. | 

The solutions show that, when disregarding a narrow layer 
adjacent to the surface of the shield, the temperature profiles 
T(t, z) in cross sections normal to the surface depend primarily 
on the thermal diffusivity k/y cp and the viscosity yu of the glass. 
This is especially true for the maximum thermal penetration 
b b(t;){mm.] since, for every one of the 35 investigated glass 
shields, both the aerodynamic heating and the melting stop about 
10 [sec.] before impact of the vehicle, and during this final 
period, k/y cp is the only relevant material parameter in the 
system of equations. 

The viscosity distribution yu(7(t,2)), which governs the melting 
process, is primarily dominated by the viscosity level, i.e., 
p because these two material parameters essen- 

T(t, s) when 
surface. The 


as is 


u* and by k/y « 


tially temperature distribution 


adjacent to the 


determine the 
disregarding a narrow layer 
thickness of the liquid film increases together with k/y cp 
seen in Fig. 2 and, therefore, the total ablation s[mm.] is more 
and more governed by k/y c, and u* as k/y cp rises. Particularly 
because melting stops about 40 [sec.] before impact, s|mm.] is 
considerably smaller than b[{mm.] as is seen in Fig. 2. The 
dependency of s{mm.] on material properties other than k/y c, 
and yw* thus has very little effect on the sum (s + }b){[mm.] 

Fig. 2 shows that b[{mm.] increases together with w*. Any 
increase of the viscosity level evidently causes a decline of 
s{mm.]. This explains qualitatively why (s + 6)[mm is 
independent of u* as seen in Fig. 2. 

It can be concluded that the thermal diffusivity k/y cp is the 
dominant material parameter for the necessary minimum 
thickness (s + b){mm.] of a glass shield, the heating and melting 
of which stops at a sufficient time before impact of the re-entry 


vehicle occurs. 
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Comment on ‘‘A Simplified Expression for the 
Period of Nonlinear Oscillations of Curved and 
Flat Panels’’ 


Robert E. Kelly 

Fulbright Scholar, Department of Aeronautics 
Imperial College, London, England 

February 24, 1960 


( _REENSPON’ has presented a simplified expression for the 
period of a panel undergoing large deformations, when the 
equation of motion may be written (using Greenspon’s notation) 
as 
wit law + Bw?) = UV 


His result comes directly from expanding the exact, elliptic 
integral solution for the period into an infinite series and integrat- 
It might be pointed out that the well-estab- 
lished first approximation of Kryloff and Bogoliuboff? gives 
nearly identical results and has much more general usage. This 
theory states that the square of the frequency for the slightly 


iag term-by-term. 


nonlinear, conservative system, 


w+ Fw) 0 


on + 


27/4) = 


1 PC 2n 
F(A sin @) sin ¢ dd 
TA 0 


where A is the amplitude and is constant for a conservative 


system. For the present case, we have 


2n 
07(A) = (a/rA fi [.4 sin? @ + (8/a)A? sin* o]do 
= a[l + (38A?2/4a)} 
The period is thus 


] 
T =2n 
Von + (38A?/4a) 


The ratio of the nonlinear to linear period is 
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TABLE | 

1 i a 5 gt be 

A/h (Greenspon ) (Kryloff-Bogoliuboff ) 
Oo! LO.) — ol 

0.5 0.891 0.889 
1.0 0.698 0.696 
1.5 0.548 0.542 
2.0 0.435 0.436 


T/T. = 
V 1+ (3BA 2, ta) 


or, using Greenspon’s formulas for 6 and a, 


a 1 
T/T. = 4! 
V: + (3/4)(6/y)(A/h)? 


where h is the plate thickness and 6 and y are constants. For 
Greenspon’s numerical example, where 6/y = 1.42, the results 
gained by use of the Kryloff-Bogoliuboff approximation are 
compared in Table 1 to those gained by use of Greenspon’s 
expansion, limiting the series to the sixth power of his defined 


RB t2., 


an | 1 
T/T. = 4 x 
Vi + (1/2)(6/y7)(A/h)? 
k? ( 9k! 225k° 
1 — + ~ 
4 64 2304 
Assuming the values of y and 6 to be typical, the first approxi- 

mation of Kryloff and Bogoliuboff certainly seems to achieve 
sufficient accuracy. This approach is, of course, much more 
general because it is valid for any form of the nonlinear restoring 


force which is dependent only on the displacement and which 
It does not first require a closed- 


’ 


is only “‘slightly’’ nonlinear. 
form solution and thus seems applicable to cases of plates which 
demand representations of the restoring force for which closed- 


form solutions are unobtainable. 
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Use of the Roughness Criterion to Refute 
Roughness as the Cause of Reported 
Transition Reversal 


Jerold M. Bidwell 

Design Engineer, Aerodynamic Applied Research, 
The Martin Company, Denver, Colo. 

February 25, 1960 


| eae lee EVIDENCE has indicated that the transition 
Reynolds number at first increases and then rapidly de- 
creases as the convective heat-transfer rate to the wall is in- 
creased (low temperature ratios, 7/7}. Although this transi- 
tion reversal has been attributed to roughness, at least one investi- 
gator! has found this effect even with surfaces smooth to 12 
microinches. It has been said? that cooling decreases the 
boundary-layer thickness and that the destabilizing effect of 
roughness becomes large in magnitude compared to the stabilizing 
effect of wall cooling. 

Braslow and Knox* have shown that the laminar boundary 
layer will be excited to transition if the Reynolds number based 
on conditions at the crest of the spherical roughness, Rx, is about 
600. The laminar boundary layer will recover and appear to be 
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unaffected by roughness if the roughness Reynolds number is 
less than about 600. This effect has been verified experimentally 
at Mach numbers up to 4 and with moderate wall cooling. 

If the criterion of Braslow and Knox yields satisfactory re- 
sults for artificially tripping the boundary layer with spheres, 
it might also be useful for estimating the degree of surface rough- 
ness responsible for the transition reversal. The roughness 
indicated by the criterion would be conservative because 
scratches would not be as effective as spheres in tripping the 
boundary layer. 

Braslow and Knox? give curves of 7 vs. R;/R, for insulated flat 
plates and cones at Mach numbers from 0 to 5. The quantity 
n is equal to (y./2x) V R, where 4; is the height of the spherical 
roughness element that is just sufficient to cause transition, x 
is the streamwise run of the laminar boundary layer, and R, 
is the Reynolds number based on conditions at the outer edge of 
the boundary layer and x. The velocity, density, and viscosity 
of the flow at the position y, in the laminar boundary layer are 
determined by the use of the laminar-boundary-layer solution 
of Chapman and Rubesin.! 

Charts similar to those of Braslow and Knox for the insulated 
flat plate* have been prepared for temperature ratios 7),/T 
from 0.08 to 1.50. To illustrate the effect of roughness height 
on laminar boundary layers, these charts have been used to pre- 
pare two examples, one with a Reynolds number, R:, equal to 
10°, and the other with R, equal to 107 (Fig. 1). The examples 
show how the ratio y;,/x varies with wall-temperature ratio, 
T/T, and Mach number for the two Reynolds numbers. The 
variation of y,./x with R, is shown for the case of 7,/7., = 1.50 
and M,, = 5.0 (Fig. 2). 

To examine this roughness effect, consider the possibility 
of transition 5 inches back from the leading edge of a flat plate 
at Mach number 2.0 where the local Reynolds number, R,, is 
taken to be 10° The insulated wall-temperature ratio for a 
laminar boundary layer at this Mach number is 1.68 and the 
wall is assumed to have been cooled to a temperature ratio of 
1.00. For this case y;/x is 1.6 X 10~* (Fig. 1), meaning that the 
spherical roughness height that will just excite transition 5 inches 
back from the leading edge is about 0.008 in. If K, was taken as 
107 instead of 10°, a spherical-roughness height of 0.0014 in. 
would be sufficient to excite transition. In most experiments of 
this kind, the surface is polished two orders of magnitude smoother 
than these values. Examination of the difference, between the 
three-dimensional 7 and the two-dimensional 7 from reference 


3, confirms the expectation that the actual difference in spherical 
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roughness height between flat plate and conical flow is relatively 
small. The fact that surface roughness is frequently of the nature 
of scratches, and that a much larger roughness of this type than 
of the spherical type is required to excite transition, is much 
more important. 

In summary, if the roughness criteria of Braslow and Knox 
is approximately valid as the wall temperature becomes very 
low, the transition reversal should not, according to these results, 


be attributed to surface roughness 
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SYMBOLS 

( = local skin-friction coefficient with mass transfer 

CH = local heat-transfer Stanton number in presence of mass 
transfer, 

ky(OT/ dy) » 
Petlelp (Ty — Ty) 

Cyo,C Ho = local skin-friction coefficient and local Stanton number 
respectively, for solid wall exposed to same free-stream 
conditions and held at the same temperature as the mass 
transfer cooled wall 


C; = specific heat at constant pressure 


’ = recovery factor 

Re, = Reynolds number = 1“,x/v, 

7 = temperature 

u, 2 = velocities, parallel and normal to surface, respectively 

ae = coordinates, parallel and normal to surface, respectively 
p = density 

v = kinematic viscosity 


* Also, Professor of Mechanical Engineering, University of Minnesota, 
Minneapolis, Minn 
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Subscripts 


0 = refers to solid wall subject to the same free-stream conditions 
and held at same temperature as the actual wall 

é = refers to conditions at outer edge of boundary layer 

u = refers to wall conditions 

r = recovery condition, refers to condition where k,y(OT/ dy) y 0 


ik IS WELL KNOWN that the magnitude of the skin friction and 

heat transfer for a flat plate in turbulent flow can be ap- 
preciably reduced if mass is released at the wall surface. Such 
a mass-transfer cooling process can be realized by the use of a 
porous surface through which the coolant is forced or by actual 
evaporation or sublimation of the surface in the presence of a 
high-temperature boundary layer. The prediction of the heat 
transfer and skin friction in such turbulent binary boundary 
lavers is of obvious technical importance. To this end, a variety 
of analytical approaches have been suggested ranging from the 
simple film theory to more complex mixing-length analyses 
The simple film theory or Couette-flow model conceptually re- 
places the actual boundary layer by a film of constant thick- 
ness, neglecting any variations in the flow direction as com 
Although this approach does 


pared with those across the film." ? 
not yield all details of the flow, it does predict the overall effect 
of mass transfer on skin friction and heat transfer quite well. 
More detailed mixing-length analyses which require empirical 
data to accomplish a final solution have been advanced for air 
injection by Dorrance and Dore,* Rubesin,* and van Driest 

The first two references are applicable to the flat-plate geometry, 
with Dorrance and Dore restricting their analysis to a Prandtl 
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number of unity and assuming the turbulent boundary layer 
extends down to the surface. The analysis of Rubesin applies 
to a Prandtl number of 0.7, considers viscous dissipation and 
allows for a laminar sublayer and a turbulent surlayer. The 
incompressible analysis of van Driest considers the laminar 
sublayer and turbulent surlayer using somewhat different as- 
sumptions than those of Rubesin for the flat-plate geometry, 
but, in addition, van Driest extends the analysis to apply to air 
injection into a turbulent boundary layer in the region of either a 
two-dimensional or three-dimensional stagnation point. The 
results of Rubesin and van Driest are essentially in agreement 
with each other for the flat-plate geometry, and their predicted 
values of the skin friction and heat transfer are somewhat higher 
than the prediction of Dorrance and Dore.* 

A number of experimental results are reported for air injection 
into a turbulent boundary layer on zero-pressure-gradient sur- 
faces,” 7 ® % 14 11 thereby providing a check on the above pre- 
dictions. Fig. 1 compares the available heat-transfer results 
in the form of Stanton numbers with the several analyses. It 
may be seen that the three analyses do not differ too markedly, 
although Dorrance and Dore are somewhat lower. The experi- 
mental data covering a range of Mach numbers from 0 to 2.7 
show considerable scatter but nevertheless the trend is correctly 
predicted by all three analyses. 

The calculation of heat transfer requires the knowledge of the 
recovery factor in addition to the above-cited Stanton number. 
The theory of Rubesin predicts that there is no effect of air in- 
jection on the recovery factor for turbulent flow. However, 
this is in disagreement with experimental evidence.” * ° Since 
some of the assumptions used in the analysis of Rubesin have not 
been confirmed, greater weight must be given to the experimental 
resuits. These are summarized on Fig. 1, and the curve drawn 
through the data may be used as a guide in predicting the re- 
covery factor for air injection into a turbulent boundary layer. 

Fig. 2 reveals the predicted dimensionless skin-friction coeffi- 
cient for the injection of several different gases into the turbulent 
boundary layer on a flat plate. Rather extensive experimental 
data for air injection” agree with the predictions of Rubesin aad 
van Driest.“\> However, the data of Mickley and Davis! show 
skin-friction vaiues much lower than the other investigators — 
even lower than predicted by the film theory. Since the majority 
of the experimental data favors Rubesin and van Driest, it 
appears at the present time that this represents the more realistic 
and certainly the more conservative estimate of skin friction 
and heat transfer in the presence of air injection. 

For foreign-gas addition to a turbulent boundary layer on a 
flat plate, the skin-friction results of Pappas!” (for hydrogen, 
helium, Freon-12), and, Ward and Harmon!’ (for Teflon) are 
available and are also represented on Fig. 2. These predictions, 


* Dorrance® discusses this discrepancy which arises from the form he 


ssumes for his velocity distribution 
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valid for low Mach numbers, demonstrate the advantage of the 


light gases, helium and hydrogen, in reducing skin friction. The 
heavier gases, Freon and Teflon, are obviously less effective in re 
ducing skin friction but surprisingly their performance is not 
appreciably poorer thaa air 

For heat traasfer, the mixing-length analysis of Rubesin and 
Pappas!‘ predicts the performance shown on Fig. 3 for light-gas 
injection into a turbulent air boundary layer on a flat plate 
The limited available experimental data for helium-air” sug 
gest an even greater reduction than predicted by the analysis 
The predicted recovery factor shows no change from the solid 
wall value when light gases are injected into a turbulent boundary 
layer and at the present time the experimental data are insuffi 
cieut to justify any definite conclusion. 

To determine the influence of molecular weight of the coolant 
gas on the heat transfer and skin friction, the mass-transfer 
parameter was plotted against the molecular weight at a constant 
value of the reduced parameter, C;/C;, and Cy/Cu, of 0.5. The 
results are shown on Fig. 4 where it is seen that the molecular 
weight appears to play a more sensitive role for the light gases, 
while the heavier gases are less sensitive to this property. An 
estimate of the reduced skin-friction coefficient and Stanton 
numbers for other coolants may be obtained from the presenta 
tion in Fig. 4. 

In closing, it may be noted that the effect of a pressure gra 
dient on the heat transfer and skin-friction coefficient in a 
binary turbulent boundary layer is apparently minimized when 
the results are presented in the coordinates shown on Figs. 1 
and 2. The analysis of van Driest for air injection indicates 
that the same curve shown for the flat plate may also be used 
for turbulent flow in the region of two- and three-dimensional 
stagnation points. 
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Some Mass-Transfer Results 
With External-Flow Pressure Gradients? 


J. R. Baron* and P. E. Scott** 
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F” SOME TIME, investigators have considered the injection of 
foreign gases into a high-speed laminar boundary layer in 
order to evaluate the effects of both the injection process and the 
behavior of specific mixtures. It is fairly well accepted that 
benefits do accrue with respect to thermal levels and heating 
rates, and that for perfect gases a low-density coolant is prefer- 
able. The findings 
relevant to the mass-transfer process in the presence of a favor- 


present note briefly summarizes some 
able pressure gradient; clearly such is of interest for blunt-body 
geometries. A detailed description of the analysis and results 
appears in reference 1 

The conservation equations for the steady two-dimensional 


flow of a binary-mixture continuum are? 








(pu)z + (pry =O | 
p(uci, + VCi,) = (pDeai,)y (1) 
plu ur + Vv ) = —pzr + (Mw Vy)y 

pcepuls +vTy) = ups + ul(yy)? + (Rk Ty)y +e Da, — hy | 


+IAI— 


Pr and Sc are Prandtl and Schmidt numbers, f = ¥ V 2s, i.e., 
a modified stream function for which f’ = u/ue, A = (pu)/(pmde, 
( ), refers to local external-stream condition, and the pressure- 


gradient parameter is 
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Here c;( = p;/p) is the ‘“‘coolant’’ mass concentration, D is the 


binary mass-diffusion coefficient, and h; is component specific 


enthalpy. The mass-transfer process requires specification of a 

mass concentration at the outer edge of the layer and a coupling 

between injection rate and surface concentration such that the 
9» 


net transfer of component ‘‘2”’ across the surface vanishes 
Introducing similarity parameters (s, 7) in accord with 


; P (pu). 4 [ p 
S = (puu). dx, = ; dy (2) 
J 0 V 2s JO Pe 
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A barred quantity implies units of free-stream value, e.g., R= 
R/R:2. Similarity solutions follow when 6, M* and Ms assume 
constant values (e.g., 4, = 0 or ©), or 8 = constant and Pr = 
Sc = 1. In view of possible blunt-body application® solutions 
with helium and air injection were obtained for several values of 
Band M, = 0 with the following boundary conditions. 


f(0) = fw, f'(0) = 0, f'"(~) = 1 
c1'(0) = (c2Scf/A)w, a(e) = 0 (6) 
T(0) = 0.8, T(o) =1 ) 

Physical property specifications correspond to 7, = S70°R. 


(see reference 1). 

Some skin-friction and heat-transfer results for helium and 
air injection are shown in Figs. 1 and 2 for the range 0.8 > -- fy, => 
O and 6 = O, 0.5, and 1.0. The ordinates are referred to their 
no-blowing, flat-plate, values. 
gradients for a flat plate and stagnation point, the negative 
pressure gradient appreciably increases both cy and St. It is 
interesting to note the overriding effect on skin friction for the 
larger 8 when helium is injected, in contrast to the usual decre- 


Over a @ range including pressure 


ment in cy. This results from an overshoot in the velocity pro- 
file (f’ > 1) which arises due to the larger accelerations imposed 
upon a low-density gas for a given pressure gradient. An in- 
crease in either 6 or f, moves the location of the maximum ve- 
locity in the profile nearer to the surface. The peak values 
asymptotically correspond to 7 = 0.8, 0.6 for 8 = 0.5, 1.0, respec- 
tively. Correspondingly, the concentration of the injected gas 
rapidly approaches unity for such 6’s (¢,,, > 0.98 for —fw > 0.8). 

Stanton number does not exhibit the same behavior for the 8 
range considered here. However, the curves do indicate a lessen- 
ing of the difference between the air and helium results with in- 
creasing negative pressure gradient. 

Boundary-layer separation is implied by cr — 0. 
asymptotic approach of ¢,, to unity with increasing /,,, it proves 


Due to the 


relatively laborious to accurately define the separation injection 
rate for 8B > 0. Nevertheless, it is clear that a substantial in- 
crease in mass flow over that for a flat plate is necessary both 
to separate the flow and to obtain similar changes in cy and St. 
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The Effect of a Deceleration Force on a Melting 
Boundary Layer 


Simon Ostrach,* Arthur W. Goldstein,** and Jesse Hamman ** 
Lewis Research Center, NASA, Cleveland, Ohio 
March 14, 1960 


Me ANALYTICAL STUDIES of ablating layers which melt be- 

fore they vaporize have been restricted to the stagnation 
region of bodies (e.g., see Sutton!). To determine conditions all 
along the body, it is important to include the deceleration force 


that opposes the downstream flow of liquid. Some treatment 


has already been given to this problem,” * but it is rather limited 
because of the idealization of the considered configuration or the 
omission of some important features of the problem. This note 
gives a more general discussion of the phenomena that occur when 
a body force acts on a melting boundary layer. The basis for 
these remarks is an analysis of the flow of a material like Pyrex 
over a two-dimensional body. The liquid density, specific heat, 

* Chief, Fluid Physics Branch. 

** Aeronautical Research Engineer 
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and conductivity are assumed constant. To determine the 
conditions under which deceleration effects are important and to 
simplify the basic equations, all variables except velocity and 
time are made dimensionless in the usual way. The character- 
istic liquid layer velocity l’* is found to be P,,h?/ Ry; by requiring 
the pressure and shear forces to be of equal order; P,, denotes 
stagnation pressure, / is liquid layer thickness, R is the body 
radius of curvature at the nose, and y; is the viscosity at the gas- 
liquid interface. Since in its initial development the liquid layer 
grows at such a rate that velocities change as a result of time- 
wise viscosity or temperature variations, the appropriate time- 
scale is found by requiring the unsteady convection and the 
necessary conduction to be of like order. The time-scale is thus 
pch?/k, where p is the density, c the specific heat, and k the thermal 
conductivity of the liquid; this time is that required for a tem- 
perature change to be transmitted through a layer of thickness 
h. Assuming now that h/R< 1 and noting that for liquids the 
Reynolds number Re = O(1), the Prandtl number Pr > 1, and 
PrU*?/cT; < 1, where 7; is the interface temperature, we can 
write the governing nondimensional equations in a coordinate 
system fixed in the interface! as 


(Ou/Ox) + (dv/dy) = 0 (1) 
3) ( Ou dp pRF, : 
= om & a TP oct ee P/ ap (2) 
oy “ oy dx Ps ; lro (x) ] I(x) 
3) h\? o7 07 
+ Pr Re V—_(x) = : (3) 
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where 7) is the distance from the body axis to the interface, t de- 
notes dimensionless time, and F, deceleration force. The simpli- 
fied form of the energy equation (3) is obtained by approximating 
the convection terms by their limiting form, v,,(07/dy), for large 
y where convection is important; near the interface, conduction 
is assumed to be dominant. 

The essential unsteady character of the problem can be noted 
from the energy equation; without the unsteady effects the tem- 
perature could not be bounded at infinity if liquid accumulated 
(v. > O) at some position. Careful examination of the pres- 
sure, shear, and deceleration forces for bodies whose surfaces be- 
come parallel to the flight direction shows that after some posi 
tion on the body the deceleration force becomes uniform and 
dominates all the other forces. The liquid layer everywhere 
downstream of this position behaves like a heated semi-infinite 
slab that is characterized by unsteady effects. These considera- 
tions mark the first important deviation from existing work on 
this problem 

Two parameters characterize the problem. The one in the 
last term of the momentum equation (2) is essentially the inverse 
of a Froude number and represents the ratio of deceleration and 
pressure (or viscous) forces; deceleration effects become signifi- 
cant when this parameter is not small. The other parameter, 
which appears in the energy equation (3), represents the char- 
acteristics of the particular liquid and gives the relative influence 
of convection and conduction. 

Liquid layer outer-boundary conditions are obtained from the 
gas boundary-layer characteristics of Cohen and Reshotko.* 
The body shape |see Fig. 3(a)] is thus limited to the class for 
which that analysis pertains. The correct interface temperature 
T; is selected to give a consistent energy balance between the 
gas and liquid layers; the velocity components and shears are 
also matched. Boundary conditions for the liquid layer thus 
are: v = —%, O0u/Oy = (R/h)(7;/Pm) at the interface (y = 0), 
and 7 = 0, « = Oat the solid surface (y = ©), where 7; denotes 
the shear at the interface. Initial conditions specified are u = 
v= T =O0att=0. Inthe calculated example, a first approxi- 
mation to 7; is used; viz., 7; = const. Although not investi- 
gated, we assume that improvements by iteration are possible 
and ignore this herein so that the other effects can be most simply 
shown. 

Taking the viscosity-temperature relation to be uy 
a; exp [b(x)y], we can write the solutions explicitly as 


u = —(e—/b)[fy + (f/b) + Fi] (4) 


/ 


r l . Y+2Z Y-Z r 
— e” erfc + erfc (5) 
Tj 2 2 2 


where 


7, = (R/h)(7;i/P on), 


Y = y/Vt, Z = Pr Re(h/R)v,, Vt 


n = Pr Re(h/R)*v,), 


At each position x along the body an iteration procedure based 
on matching the assumed viscosity-temperature law with the 
actual one for Pyrex! at two points is used to determine b(x) and 
Gatti 

The liquid layer development with no deceleration can be seen 
in Fig. 1. Negative values of v~ indicate a thinning of the liquid 
layer. Thus, liquid flows away from the stagnation region and 
tends to accumulate downstream because the shear and pressure 
decrease. At some later time the thicker layer results in a 
greater pressure force, and the liquid is moved rearward in a 
wavelike pattern. In proceeding downstream from the stagna 
tion region with a deceleration force, however, we found that 
after some time our calculation procedure failed at some point 
It was felt that, among other reasons, this occurred because the 
forward integration could in no way be influenced by upstream 
flow due to the deceleration force. Accordingly, asymptotic 
(in x) forms of the basic equations were integrated by the same 
procedure starting at the rear of the body. The effect of the 
deceleration force (for pF,R/P» = —0.1) can be seen in Fig. 2. 
The liquid no longer continues downstream but accumulates at 
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some point on the body. The inhibition of downstream liquid 
flow results in more effective heat shielding by the layer. It is 
noteworthy that fortuitously v, computed from the front and 
the rear for 1.45 sec. match smoothly. Representative profiles 
at different body stations for this case are shown in Fig. 3. Flow 
reversal and upstream flow at the back end due to the body force 
are clearly evident from the velocity distributions. Although 
the temperature profiles look quite similar at this relatively short 
time, their slopes at the interface indicate relatively higher heat 
transfer into the layer in the stagnation region ( = 0) with the 
lowest heat transfer at the position of maximum v,(& = 0.6) 
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SYMBOLS 


E = elastic modulus 
E, = reduced modulus 
E; = tangent modulus 
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I = generalized inelastic modulus en a — 
I = moment of inertia TABLE | 
. = length ae ae 7 a ” 
: iedee Column Type M, M; Ger 
M, = external bending moment SS 
M; = internal bending resistance perfect elastic Pw — EI (d2w/dx2) rE p? L? 
P = axial compressive load imperfect elastic P(w+w) —EI(d*w/dx?) w?Ep?/L? 
w = lateral deflection : perfect inelastic >w —E,I(d?w/dx?)  2*E,p?/L? 
= initial lateral imperfection imperfect inelastic P(w+w) —Ed(d*w/dx?) x2*E:p?/L? 
p = radius of gyration ee i“ 
qa = axial compressive stress 
ocr = compressive buckling stress . i ” F : 

For perfect columns, #@ = 0, and therefore the solution is deter 


INTRODUCTION 


Re gaeapige RESULTS on aluminum-alloy columns led 
Shanley! to re-examine the basic assumptions of inelastic 
column theory in a now classic paper published in 1947. He 
concluded that if axial straining and bending proceeded simul- 
taneously at buckling, no strain reversal would occur and, there- 
fore, the tangent modulus is the correct effective modulus for in- 
elastic buckling of a perfect column. Shanley’s analysis of an 
idealized type of inelastic column indicated that a continuous 
spectrum of possible bent configurations exists for which the 
lateral deflections increase from zero at the tangent-modulus 
stress to infinity at the reduced-modulus stress. Thus the tan- 
gent-modulus stress is the lowest value at which a bent configura- 
tion becomes stable and, therefore, can be considered as the 
buckling stress. 

In commenting on Shanley’s contribution, von Karman! indi- 
cated that the reduced modulus is correct when the stability 
analysis is based on the requirement that the axial load remain 
fixed at the exchange of equilibrium configurations. He con- 
tended that because of the nonreversible character of inelasticity, 
it is necessary to revise the definition of the stability limit. Thus, 
for inelastic buckling, we should seek“ . . . the smallest value of 
the axial load at which bifurcation of the equilibrium positions 
can occur, regardless of whether or not the transition to the bent 
position requires an increase of the axial load.’’ From this 
definition, the tangent-modulus buckling stress originally pro- 
posed by Engesser is the correct solution. 

The analysis presented in this paper utilizes Shanley’s con- 
cept that axial straining and bending proceed simultaneously 
without strain reversal in the development of tangent modulus 
theory. However, it is not assumed that this concept applies to 
a perfect inelastic column, nor is it required that the definition of 
the stability limit be revised for inelastic buckling. Instead, the 
theory of an imferfect inelastic column (imperfect in a geometric 
sense) is developed in which case the assumption that axial 
straining and bending proceed simultaneously without strain re- 
versal is completely rational. 


ANALYSIS OF COLUMN TYPES 


Since the pertinent differential equations involve equilibrium 
of bending moments, we shall be particularly interested in the 
external bending moment, J/,., and the internal bending resist- 
ance, /; for the four column types: perfect and imperfect elastic 
columns, perfect and imperfect inelastic columns. Only in- 
finitesimal lateral deflections are considered because we are con- 
cerned here with the buckling behavior. 

The usual relations for the bending moments are listed in 
Table 1 where w is the lateral deflection of the column always 
measured from the initial position, #. For perfect columns, of 
course, w = 0. 

Although the values of /; given in Table 1 are familiar, it 
should be noted that for the perfect inelastic column, E; is based 
on the assumption that bending starts at buckling under a fixed 
axial load and, therefore, a strain reversal occurs. On the other 
hand, for the imferfect inelastic column, axial loading and bend- 
ing proceed simultaneously from the start due to the geometric 
imperfection. Hence, the tangent modulus is appropriate in this 
case. 

In all cases, we seek a solution of the equilibrium equation 


(d*w/dx?) + (o/Ep?)(w + w) = 0 (1) 


mined for characteristic values of o as given in Table 1. For 
imperfect columns, the initial geometric imperfection can be con 
veniently represented by 


© 


Similarly, the lateral deflection for a simply supported column 
w= z. w, sin (n x x/L) (3 
n=1 
By carrying out the analysis for the imperfect inelastic column, 
it can be determined that 
" (n?27?E.p?/L?) 
W,/tin = = (n Reee<c.) £6 
Oa 
As the axial compressive stress o, approaches 7?E,p?/L?, which 
can be identified as the limiting value of axial compressive stress, 
ocr, for this case, the m = 1 term dominates Eq. (3). In addition 
as oq approaches o-,, the tangent modulus becomes substantially 
constant. 
Thus, in the region of buckling, the lateral deflection at the 
center of the column is approximated by 


w/w, = [(aer/oa) — 1]7! (5) 
where Cer = w*E.p?/L? (6) 
For the imperfect elastic column, EZ, = £ and Eq. (6) reduces to 


the familiar Euler equation as indicated in Table 1. 


COMPARISON OF SOLUTIONS 

At this point, it is most interesting to compare the solutions for 
perfect and imperfect columns. In the elastic case, the limiting 
axial compressive stress for a column containing small initial 
geometrical imperfections is identical with the buckling stress 
of a corresponding perfect column as shown in Table 1 

In comparing the inelastic solutions given in Table 1, however, 
the buckling stress of a perfect column is always higher than the 
limiting axial compressive stress of an imperfect column, Eq. (6) 
Since the latter represents the buckling stress of a column with 
vanishingly small initial imperfections, the tangent-modulus 
buckling stress can be considered as the correct theoretical in- 
elastic generalization of the Euler stress. On the basis of experi- 
mental evidence, it has so been considered for a long time. 
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Reynolds-Analogy Parameter for the Laminar 
Boundary Layer With Blowing, Pr = 3/4 


C. R. Faulders 
Aero-Space Laboratories, Missile Division, 
North American Aviation, Inc., Downey, Calif. 


March 22, 1960 


SYMBOLS 


h = enthalpy 
k = thermal conductivity 
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Fic. 1. Influence of blowing on Reynolds-analogy factor. 
Pr = Prandtl number 
7 = temperature 
u = x-component of velocity 
v = y-component of velocity 
x = coordinate parallel to surface 
y = coordinate perpendicular to surface 
uw = coefficient of viscosity 
p = mass density 


p | ‘HE SHEAR STRESS and the heat-transfer rate at the wall are 
related to the skin-friction coefficient, Cr, and the Stanton 


(> Cr a 1 
Hu = Pata° (1) 
Ov) u 2 


oT ’ 
ky, wef J = Cypot.(haw — hw) (2) 


The subscripts w and © denote wall and free-stream values, 
respectively, and haw is the adiabatic wall enthalpy. With a 
Prandtl number of unity Cy is exactly equal to Cr/2, while for 
Prandtl numbers other than unity and a laminar boundary 
layer without blowing or suction the ratio 2Cq/Cr can be closely 
In view of the import- 


number, Cy, by 


approximated by (Pr)~?’, reference 1. 
ance of mass injection into the boundary layer as a means for 
reducing heat transfer, the variation of 2Cy/Cr with blowing 
rate is of interest. The analysis presented here is restricted to 
the injection of air into air. 

The flat-plate laminar boundary equations 


(pu )z + (pv)y = O (3) 
puu, + ptuy = (puy,), (4) 
puh, + pth, = (RTy)y + wu,’ (5) 


can be reduced to ordinary differential equations following the 
usual procedure of introducing the stream function 


}UpQXp 4s . 
y = y oe i (6) 
Pe 
and the similarity parameter 
1 lu , F p 
n=- ora Sd ds (7) 
2 Xe 0 Po 


For (pu) constant and equal to (p,.u..), Eqs. (6) and (7) permit 
the reduction of Eq. (4) to the Blasius equation, 
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Primes denote differentiation with respect to 7. Enthalpy can 


be substituted for temperature in Eq. (5) using 


dh c,pdT (9) 


F 


where c, is the specific heat, as long as enthalpy is a unique 


p 
function of temperature 


a function of u only, introducing the Prandt! number, assumed 


Considering now the enthalpy to be 


constant, and using Eqs. (4), (6), and (7), the energy equation 
becomes 
(d*h/du*) + ((1 — Pr)0(dh/du)| + Pr 0 
Au) — ((2/u f/f" (10) 


The coupling between the enthalpy and velocity distributions 
when Pr ~ 1 is represented by the function 6. With boundary 


conditions 


u ten: Aah 


the first integral of Eq. (10) evaluated at the wall is 


dh (h, — h, 
j —79e + ou, (11) 
du) 


mit —_ tm Ug. 
| e~ Mo 95 dy 
J0 


where o represents other terms that depend only on the Prandtl 
number and the velocity profile. Recognizing that h, Risa 
when the derivative of enthalpy at the wall is zero, the right- 
hand side of Eq. (11) is reduced to just the first term, but with 
(ha — hw) replaced by (haw — hw). Combining Eqs. (1), (2), 


(9) and (11), 


Cu l ] 
ae 9 ~ P eu ( 12) 
F/< 1 > G 
| e~( Pr So “0ds ay 


Eq. (12) has been evaluated for P; 3/4 and values of fy, 
representing the blowing rate, varying from zero to —1.2. The 
function 6(2) was obtained from tabulated solutions of Eq. (8).? 
The numerical results are plotted in Fig. 1. It can be shown 
that 2Cy/Cr approaches infinity as f, approaches a value of 
— 1.238, corresponding to incipient separation, as long as the 


Prandtl number is less than unity 
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Rarefaction Effects in Low-Speed Turbulence 
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SYMBOLS 


d = diameter of hot wire 
t'/F = root-mean-square fractional fluctuation in hot-wire heat- 
transfer rate 
zg = Olog Nuq/2 log (1 + 7) 
M = Mach number 
* Professor of Chemical Engineering 
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TURBULENCE KNUDSEN NUMBER A 2 


Fic. 1. Correlation of turbulence intensities by a turbulence 


Knudsen number. 


n* = Olog Nu~*/O0 log Rey 

Nu = Nusselt number (hd/k) 

Rep = Reynolds number of turbulence generator 
u’/U = relative turbulence intensity 

» Kolmogoroff scale (v3/e)!/4 

€ rate of dissipation of turbulence energy 
rd molecular mean free path 

P = kinematic viscosity 

T = temperature loading (A7/7) 

ta) 0 log ¥/0 log Nu 

y Nuow/Nu, wire-length correction 

( da refers to an infinite-length wire 

( )* = refers to continuum condition 


— MEASUREMENTS were made at constant Reynolds 
number of the turbulence intensities at the center of the 
flameholder plane of a 6-in. model ram-jet in cold flow of air at 
pressures varying between 1.0 and 0.04 atm. Three levels of 
Reynolds number based on ram-jet diameter were studied 
20,000, 35,000, and 57,000. The maximum velocity employed 
was 273 ft./sec. As the test pressures were decreased below 
0.2 atm., the turbulence intensity at constant Reynolds number 
decreased significantly. 

The following equation was developed and used to compute 
the response to turbulent velocity fluctuations of a nonideal hot 
wire operating in isentropic, low-speed flow: 

"0 = (Lo + 7g). f'/F (1) 
' 1 — AyNu(r/d) n* 
An equation suggested by Collins and Williams! was used to 
correct for the slip-flow effect on the response of the hot wire: 


1/Nu,, = 1/Nu,* + Ad/d (2) 


The second term is proportional to the temperature jump. The 
coefficient A was determined from mean heat-transfer measure- 
ments taken simultaneously with the turbulence measurements. 

The turbulence-intensity results calculated by means of Eq. 
(1) from measured values of f’/F are presented in Fig. 1 as a 
function of the turbulence Knudsen number \/n. This method 
of treating the results eliminates a separate effect of Reynolds 
number. The length 7 is the Kolmogoroff scale characteristic 
of the eddies responsible for the viscous dissipation of turbulent 
energy, and may be considered a measure of the lower bound of 
eddy sizes existing in the system. The correlation of intensities 
over a range of Reynolds numbers and pressures by a single 
function of the turbulence Knudsen number strongly indicates 
that the rarefaction effect is a result of the noncontinuum nature 
of the fluid. 

These results are a refutation of a fundamental premise of tur- 
bulence theory: that for all practical conditions the fluid may 
be regarded as a continuum, and hence the validity of the Navier- 
Stokes equations may be assumed, even when applied to the 
smallest scale components of the motion. This assumption is 
undoubtedly correct at normal atmospheric conditions, but it 
has been demonstrated here that significant rarefaction effects 
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do occur in practical systems at pressures far higher than pre- 
viously have been suspected. 

The initial dissipation turbulence Knudsen number is ex- 
pressible in terms of the Mach number and a Reynolds number, 
Rep, based on a characteristic dimension of the turbulence- 
generating device 

A/n = AM/(Rep)'4 (3) 


The coefficient A is presumably a function of the drag coefficient 
of the generating device and of the system geometry. Ina flow 
regime where the drag coefficient is independent of Rep, the ratio 
M/(Rep)'!* may replace \/n, and this is the parameter which 
determines the existence and magnitude of rarefaction effects on 
turbulence dissipation rates in geometrically similar systems 
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the Eccentricity and the Perigee Angle 
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SUMMARY 
A graphical method is presented for the determination of the 
eccentricity and perigee angle in the two-body problem for a 


given initial velocity and angle. 


SYMBOLS 


e€é = eccentricity 

v = velocity 

V = dimensionless velocity, v/v 

a = radial component of velocity 

8 = horizontal component of velocity 

y = velocity angle with respect to horizontal 
6 = perigee angle 


Subscripts 
c = circular 
0 = initial condition 


_— ECCENTRICITY AND PERIGEE ANGLE can be derived from 
the usual simplified two body assumptions as 
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Fic. 1. Illustration of graphical method. 
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Fic. 2. Eccentricity and perigee angle. 


e= [a?B? + (sp? — 1)?]’ 
@ = cos 1 [(p2 — 1)/e] (1) 
where 
a = Vosin yo 
B = Vocos yo 


The perigee angle as given in Eq. (1) is independent of whether 
the velocity is applied in an upward or downward direction. 
To avoid ambiguity, the perigee angle will be considered as 
measured in the opposite direction from the rotational direction. 
Downward applications of velocity will thereby result in perigee 
angles between 180 degrees and 360 degrees, due to the initial 
point being beyond the apogee point. 

Fig. 1 shows several initial velocity vectors represented on 
an a — 8 graph. A line can be drawn from where the initial 
velocity vector intercepts the point 6? plotted on the 8 axis to 

The horizontal distance between the 
The vertical displacement of the inter- 


the point a = 0, 8B = 1. 
two points is |1 — 6? 
section is found from symmetrical triangles to be a. 

It is thereby seen that the above constructed line represents 
the eccentricity with the perigee angle being the angle between 
the line and the horizontal axis (see Fig. 1). 

The following steps can be used in the graphical determination 
of e and @. (1) Draw the line representing the initial velocity 
vector. (2) Locate 6? on the 6 axis. (3) From where this 
point intersects the velocity vector line, draw a line (represent- 
ing the eccentricity) to the point a = 0, 8 = 1 (the perigee angle 
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is measured from the 8 > 1 portion of the abscissa in a counter- 
clockwise direction). 
The perigee angle and eccentricity are shown in Fig. 2 


SPECIAL CASES 

(l)B@=1,".e2z a 

The perigee angle is 90 degrees for upward Vo, 270 degrees 
for downward Vp. The semimajor axis is perpendicular to the 
initial radius vector. 

(2) a =0, .. e = B? — 1 

The perigee angle is zero degrees for 8 > 1, 180 degrees for 
6 < 1, corresponding to a perigee and apogee at the initial point, 


respectively. 


(3)B =0, .e =1 
A vertical orbit with a perigee angle of 180 degrees. 
(4) a? +B? = 2, .e = 1 


A parabolic orbit with perigee angle given by 
6 = cos '(B? — 1) = cos~!(1 — a?) 
(5) a? + B? = 1, .. e = (1 — B? + £B*) ”? 
Conditions for initial velocity vector parallel to the major 
axis, obtained by equating the slope of the velocity line to the 
minus tangent of the perigee angle 
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SUMMARY 

w pen STRESS FIELD, displacement field, and torsional rigidity 

which are presented for a prismatic bar of hexagonal cross- 
section have been obtained by means of the Galerkin technique. 
For purposes of mathematical analysis of the problem only 
two terms of the stress function have been used. The con- 
vergence of the solution is demonstrated, and the resulting 
stress distribution between opposite corners and between opposite 
sides is tabulated. The analysis indicates that a maximum 
shear occurs just inside the boundary from the midpoint of a flat. 


FORMULATION OF THE PROBLEM AND ITS SOLUTION 
The use of prismatic bars for torque transmission is wide- 
spread, and a knowledge of their elastic torsional behavior is 
important, especially in structural design of modern aircraft 


engines. 
To date, however, the approximate solutions for only a few 
types of sections have been obtained [1, 2, 3, 4]. The solution 


of the torsion problem for a bar with hexagonal cross-section can 
be approached through the Ritz, Kantorovich, or Galerkin 
methods. In this paper, the Galerkin method is adapted for the 
mathematical analysis of the problem 

The Prandtl stress function W is determined by finding the 
solution for Poisson’s differential equation 


L{[Wv(x, y)] =V2V (x,y) +2 = 0 (1) 


in the region S, and subjected to the boundary conditions, 
v(T) = 0, on the boundary of S. (See Fig. 1.) 


Assume 
N 
V(x,y) = 2 Amnm(x,y); (m = 1,2, ,N) (2) 
m=1 


Then the Galerkin method of stated problem implies that 


f L [W (x, y)] ar (x,y) dS = 0; (r = 1,2, ,N) (3) 
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where 
. y ow(T) =0 (5) 
A 
. and 
\ 
N Pm (x, y) = Ar + Ao(x? + y?) + Asxy + ; (6) 
\ i, — If only two terms are used, Eqs. (2), (4), and (6) lead to 
W (x, y) = Aim (x, y) + Aone (x, y) (7) 
Evidently 
x er, ae a m(x, y) = w(x, y) (8) 
no (x,y) = (x? + y?) w(x, y)) si 
Eq. (3) by virtue of Eq. (7) becomes 
Fic. 1. Torsion of a hexagonal beam. 


ff {V2[Aim(x, y) + Aon(x, y)] + 2} m(x, y)dS = 0 (9) 


and 


ff {V2[Aim(x, y) + Aane(x, y)] + 2} m(x, y) dS =O (10) 
s 


Moreover, for a hexagonal cross section with side of length r 


where nm (x,y) are coordinate functions, and A,, are arbitrary 
constants to be determined from Eq. (3). However, Eq. (2) 
can be written as 


W (x, y) = w(x, y) Pm (x, y) (4) a 
w(x, y) = [y? — (3/4)r?][y? — 3(x + 7)?] [y? — 3(x — r?)] (11) 


— . From Eqs. (9) and (10) it follows that 


aff mts ff nv*nds ~~ 2 f f nas ff nvnas 


A, = = (12) 
f forivnas ff nvnas _ ff nvinds f f mV *nod S 
of finds ff nvnas - 2f fas ff nvenas 
s Ss Ss Ss (13) 


Ao = - 
ff nvenas ff nx mas 


Hence, the values of A; and A» are obtained as soon as a set of finite integrals is evaluated. 
the surface under consideration, Fig. 1, one obtains 
3 3 
Vv‘ ff m aS Vv f ffir dS 
4 { 
V3 aie si ». V3 —_—* oe 
: mV *m dS = —8.646 X 10r!2;) _ nv *no dS = 1.686 & 10?r%6 (14) 
47. é 


V/3 V/3 : 
! 5 mV *n dS = —4.164 X 1074; mV *n2dS = —1.094 & 10r'4 
4 4 


Using Eq. (8), and after integrating over 


— 2.97478; —7.06 X 107}710 


From Eqs. (12), (13), and (14) it follows that 








A, = —6.77 X 10-*7-4; A. = —8.36 X 10-37r-4 (15) 
Hence ow(x, y) | 
ates " o a - O 0 Ga = 
W(x,y) = —6.77 K 10-2/r*; w(x, y) [1 + 0.124(x? + y?)r/?] (16) oy 
Note that x? + y® <r’, and therefore the series expression for 0 0 G Ov(x, y) (17) 
A rij = = 
the analytical structure of W(x, y) converges very fast, as to be ‘ si ox ‘ 
expected. This also shows that for many engineering problems  aw(x, y) — aW(x, y) 
using only the first term in approach described above yields a Ga - —Ga > 0 
x 


very good qualitative result. 


The stress tensor is given by where G is shear modulus and a twist per unit length. Evidently, 





TABLE 1 
The Value of Shear Stresses for Ratio y/r and x/r, Respectively 





0.1 9:2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
(O, y/r) oma os, = sO 6 AR mae 
— SZ 0.170 0.332 0.478 0.599 0.687 0.740 0.746 0.703 ; 
Gar 
Tzy(x/r, 0) i ais i = ee oe = a ss 
—— 0.172 0.333 0.478 0.595 0.675 0.708 0.672 0.557 0.341 0.0 
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ae » Ga J w(x, y) x? + y? 
T22 —60.44 << 10 - 7 1 + 0.124 - + 
rt} oy ta 
2y 
a(x, y) 
r? f 
ao , Ga fdw(x, y) x* + y* 
tT: Git xX Or — « 1 + 0.124 F a 
r§ lox r? 
2x 
a(x, vy) 
iad f 
(18) 


The displacement field is characterized [2] by 


— ays 
ui; = | axe (19) 
agp(x, ¥) 
where 
(x ,0) dy (x,y) oy 
o(x, y) = dx — dy + const. (20) 
(0,0) Ox (x,0) Ox 
with 
W(x, 9 (1/2)(x? + y?) — [(6.77 XK 10-*)/r*] w(x, y) X 
[1 + 0.124(x? + y?)/r?| (21) 


Finally, the torsional rigidity D, is given by 


D, = 2G f five. y) ds (22) 
Hence 
8G (V3 r/2) (r-—(9/V/3) 
D, 7 [Aim(x, vy) + Aone(x, y) |dxdy 
V3J0 J0 
D, = 0.962 Gr? (23) 


Therefore, 
D,, = (2/7r)Do 


where Dy is torsional rigidity for a bar with circular cross section 

Calculated values of 7. versus y/r and 7,, versus x/r are shown 
in Table 1 and are plotted in Fig. 1. A maximum shear is 
indicated for 7., at y/r ~ 0.7; physically, however, a maximum 
should occur at the boundary, y/r = 0.866 
the approximating 


Evidently this 


discrepancy is a result of nature of this 
solution, and better agreement would result if additional terms 


were used in the stress function 


~ 
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Periodic Temperature Distribution 
in a Two-Layer Composite Slab 


W. F. Campbell 
National Aeronautical Establishment, Ottawa, Ont., Canada 
March 28, 1960 


i» A RECENT CONTRIBUTION to the Readers’ Forum, under the 
above title, Stonecypher! outlined a method for finding the 
periodic temperature distribution in a two-layer composite slab, 
one exposed surface of the slab being insulated and the other 
subject to a sinusoidal temperature variation. Perfect thermal 
contact between the two layers, and constant thermal properties 
were assumed. 

Two years ago I drew attention in these pages to a method for 
determining the transient temperature in such a two-layer slab 
resulting from a triangular heat-input pulse.?:*»* I should like 
to point out that this same method also is applicable to the case 
where one external face is given a sinusoidal temperature vari- 
ation with time. The method is based on the analogy between 
one-dimensional heat flow and the flow of an electric current 
in a simple transmission line having only series resistance and 
parallel capacitance 

Briefly, a temperature ‘‘wave’’ may be considered to be propa- 
gated from the external face of the first layer toward the inter 
face with the second layer. At this interface the wave is par 
tially transmitted and partially reflected, the proportions de- 
pending on the thermal properties of the two layers. The trans- 
mitted wave is then totally reflected at the insulated face of the 
second layer, while the wave reflected back into the first layer 
is reflected again at the outer face of the first layer with a reflec- 
tion coefficient of —1. These two waves then approach the 
interface and, again, are each in turn split into two partly re- 
flected and partly transmitted waves, and soon. Rapid attenu 
ation with the distance a wave travels ensures that for many 
problems of interest the contributions of only a very few waves 
must be summed to give the temperature at any particular 
point 

The solution for a semi-infinite slab all of one material with a 
Q is, for an internal 


temperature 7 = 7) cos wf at the face x 


point x 
£ cos (wl — 


T Tye 


&) 


where &— is a nondimensional distance depending on the angular 


frequency w and the thermal properties of the material of the 


slab® 
REFERENCES 
Carter, W. J., and Oliphint, j. B., TVorston of a Circular Shaft with bs xV w/2a 
Diametrically Opposite Flat Sides, J. Appl. Mech., Vol. 19, No. 3, Septem a k pe 
ber, 1952 * es 
wi ; k thermal conductivity 
Shirely I K., and Stanisi¢, M M , On the Torsion of ( ylindrical Beam , thermal capacity 
Having a Trapezoidal Cross Section, Journal of the Aero/Space Sciences 4 
Readers’ Forum, Vol. 26, No. 7, july, 1959 p density 
3 Panov, D. Y., Concerning the Torsion of Nearly Prismatic Rods, Prikl = * . : 7 
Mat. Mekh., Akademiya Nauk SSSR, New Series, Vol. 2, pp. 159-180, 1938 rhe solution for the two-layer slab with the first layer extend- 
Doklady Akademii Nauk SSSR, New Series, Vol. 20, pp. 251-253, 1938 ing from x = Otox d,, and the second laver from a d, to 
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Readers’ Forum, Vol. 27, pp. 461-463, June, 1960 T = T, cos wt at x 0) 
OT /Ox 0 at v dj +d 
= ‘ = and temperature and heat flow continuous across x d;, may 
then be written down as 
o . - ( — tt; —2a(p+q) —2b(qa+r) 
N= T dX YX YY (Re) SS) Ray(—)? +4 Mh (p, g, New ~ 24 +O — 24” x 
p=0q=07=0 
cos [wt — & — 2a(p + q) — 2b(¢ + 7)] — Milp, qr yet f1—2a(P+9)—2b009+7) Eos [wt + & — 2a(p + q) — 2b(¢ + 1r)]} 
1 T  B } 2 (Riz)?(Si2)4 7 Sar) R»)'(—)P ta} Me(p, q,’ je & — 2a(p + a) — 2b(a +1) x 
p=0q=0r=0 
m P ; ‘ , 4-fo— 2, a+1)—2b(q+ fi - ‘ l ‘ 
cos [wt — & — 2a(p + g) — 2b(q + 1r)] + Nol(p, 9, ret? “APT 1) —20(a+r +1) eos lwt + & — 2a(p +q + 1) — 2b(¢+7r4+1)]} 
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where the distances are nondimensionalized 


/ el F ; 
a= d,V w/2a, b= dV w/2a2 
i= nV w (2a & =~at+(xm— by 03 2a2 


1 


reflection and transmission coefficients have been introduced 


Ry = (m2 — 1)/(m2 + 1) = —Ra 
Sio = Zure/(mi2 + 1) = mre So 
m2 = V kipics/Repr2ce 


and ‘‘multiplicity factors” are defined as 


M(p, q, r) = the number of different and distinct paths from 
the plane of the temperature source at x = 0, to the plane 


x, arriving at x in the direction of x increasing, which involve, 


in any possible order, p reflections back into layer No. 1 
at the interface x = d,, g passages through the interface in 
the direction of x decreasing, and r reflections back into 
layer No. 2 at the interface 

N(p, q, 7) = the number of different and distinct paths from 
the temperature source at x = 0 to the plane x, arriving 
at x in the direction of x decreasing, which involve, in any 
possible order, p reflections back into layer No. 1 at the 
interface x = dj, g passages through the interface in the 
direction of x decreasing, and ¢ reflections back into layer 
No. 2 at the interface. 

The subscript 1 refers to points in, or properties of, the first 

laver, and similarly, subscript 2 refers to the second layer. 
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A Method of Reducing Aeroelastic Effects 
of Highly Swept Wings 
Edward J. Hopkins,* Don W. Jillie,* and Raymond M. Hicks* 


National Aeronautics and Space Administration, 
Ames Research Center, Moffett Field, Calif. 


March 24, 1960 


SYMBOLS 


Che = lift-curve slope based on projected plan area, per deg. 

dC p/dC 2 = drag-rise factor 

CDo = minimum drag coefficient based on projected plan area 
L/D)max = Maximum lift-drag ratio 

dC»/dCy, = aerodynamic center location from 0.35¢ in terms of ¢ (nega 

tive values for rearward locations) 
A = aspect ratio 
b = wing span of flat wing. in. 


mean aerodynamic chord, in. 


ins 
i} 


modulus of elasticity, lb. /sq. in. 
To = moment of inertia of the flat wing at section Jo — Jo of 
Fig. 1,in.4 
M = Mach number 
q = free-stream dynamic pressure, Ib./sq. in 


ghee GREATEST EFFICIENCY for a lifting surface at supersonic 
speeds, according to the theoretical considerations of refer- 
ence 1, can be attained if the leading edge is swept well behind 


* Aeronautical Research Scientist. 
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Fic. 1. Estimated effect of flexibility on the aerodynamic 


characteristics of a highly swept arrow-wing with full leading- 
edge thrust assumed; J = 3.0. 


the Mach cone and the highest aspect ratio which is structurally 
possible is employed. Such a wing, designed for a Mach number 
of 3.0, would have 80° of sweepback.? Aeroelastic effects have 
been shown’ to be considerable for a wing with 60° of sweepback 
and designed for a Mach number of 2.0. The wing? shown in 
Fig. 1 was found theoretically to have considerable loss in maxi- 
mum lift-drag ratio attributable to aeroelasticity. This wing 
has 12-per cent-thick Clark-Y airfoils normal to the wing leading 
edge. If it were of solid aluminum and flying at a dynamic 
pressure of 2,400 Ibs./sq.ft. (flexibility parameter gb*/EJ) = 7.8), 
analysis indicates that the wing would deflect so as to reduce 
the maximum lift-drag ratio about 30 per cent 

In an attempt to improve the aerodynamic characteristics, 
calculations were made for this wing with each panel warped to 
conform with the curvature of a cylinder having its axis parallel 
to the flight direction (Fig. 1). Stiffening the wing in this manner 
would result in no change in local angle of attack and only small 
changes in minimum drag. 

The linear theory of reference 4 was used to calculate the sur- 
face loading for the flat rigid wing. Surface loadings for the flat 
flexible and the warped flexible wings were derived from those for 
the rigid wing by an iteration process involving load adjustments 
for the calculated wing deflections. No more than four approxi- 
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mations were found necessary to arrive at the final loading curves. 
In all cases the wing panels were treated as simple beams. Values 
of the lift and pitching moments were obtained by mechanical 
integration of the loading curves 
The following assumptions were made in the calculations: 
a) Modulus of elasticity is constant along the wing span 
(b) Airfoil taper of the flat wing is uniform along the span. 
(c) The warped wing is derived from the flat wing by warping 
each panel to conform with the portion of a circular cylinder 


shown in Fig. 1 

(d) A negligible deflection occurs ahead of Io-Ip in Fig. 1. 

(e) Center of load at each spanwise station is coincident with 
the neutral axis; thus zero twist about the neutral axis is as- 
sumed. (This assumption is considered valid because the analy- 
sis indicated that the maximum change in local angle of attack 
due to twisting was less than 1 per cent of the value due to 
bending. ) 

Although the curvature of the warped wing is not considered 
to be optimum, the wing warpage was so effective that no signif- 
icant wing distortion or change in the aerodynamic character- 
istics due to aeroelasticity was calculated for values of gb*/E/ 
up to 7.8 (see Fig. 1). For small values of gb*/Elh, the warped 
wing has a slightly smaller (L/D),,,, than the flexible flat wing 
because of the loss in lift associated with the smaller projected 
irea. 

It was noted in the analysis that for a given sized wing, as the 
dynamic pressure was increased, more of the wing tip of the flat 
wing became completely unloaded; a result which was unaltered 
by an angle-of-attack change within the range of validity of linear 


theory 
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Higher-Order Corrections to Taylor’s Solution 
for Weak, Normal, Shock Waves* 


Martin Sichel 

Assistant in Research, Gas Dynamics Laboratory, 
Princeton University, Princeton, 

March 28, 1960 


Te FIRST ORDER THE SOLUTION of the one-dimensional Navier- 

Stokes equations by the method of expansion in a small 
shock-strength parameter is identical to Taylor’s solution for 
weak shocks 
solution have been found from the Navier-Stokes, Burnett, and 


Though higher-order corrections to Taylor’s 
Grad Thirteen-Moment Equations, !: ?: * these solutions either 
are slowly convergent,! or fail to give the shock structure directly 
as a function of displacement x, and the fluid properties.” 
Hence, rapidly convergent second- and third-order corrections 
are obtained below by direct substitution of series for u, p, 7, and 
p into the Navier-Stokes or continuum equations‘ which, accord- 
ing to recent shock-structure measurements, should be valid for 
M, <2. 

Dimensionless velocity, u; temperature, 7; longitudinal vis- 
cosity, #”; the shock-strength parameter, e; and reference quanti- 
ties are defined below. Subscripts 0, 1, and 2 refer to reference, 
upstream and downstream quantities, respectively, and dimen- 

* This research was sponsored by the Office of Scientific Research, USAF, 
\ir Research and Development Command, under Contract AF 49(638)-465 
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FORUM 6: 


sional quantities are indicated by the subscript d. Hence, there 


follow: 
Uoa = (1/2)(tia + Mea); & = Ua/tod; w= pa” /p 

€ = (ta — Uea)/(tia + Wea); poa” ba’ (Toa); k ka/k 
Toa = (toa?/yR)[1 — &(y + 1)/2); 
0 = TayR/uoa? + &(y + 1)/2 


After the pressure and density have been eliminated with the 
perfect-gas and continuity equations, the momentum and 


energy equations become 


O/y + u? — eu” u(du/dx) — &(y + 1)/25 u(y + 1)/3 (1 
—60/y(y¥ — 1) + u?/2 + [k/Pr’(y — 1)]e(d0/dx) — 
e(y + .1)/27 = uly + 1)/7 — (y¥ + :1)/2%y - 1) (2 


Now uw, 6, and w” are expanded in the series 


uw = 1 t+ eu(x) + c2u(x) + cu rv) + 3 
0 1+ @)(x) + 6) (x) + Q(x) 4 j 
a 1 + eu Dix) + e*u (x T €p x T 5 


Since Pr” is assumed constant, there follows: 
a” k 6) 


The coefficients, wu“, of the series for xn” can be expressed in terms 
of the temperature coefficients 6° if a viscosity-temperature law 
of the form pa” ~ Ta is assumed 

The expansions are substituted in the conservation equations 
(1) and (2 Equating coefficients of like powers of ¢€ vields 
simultaneous equations for 6 and wu from which either @ 
or u\ may be eliminated. The nonlinear equation for u\ is 
identical with the steady-state Burger’s equation® with the 


solution 





uo — tanh Ax r 
a) = (y — 1) tanh Ax 8 
1 ( + 1)/2]/{1 + — 1)/P 9 
TABLE 1 
€ 0.10 an 0.5 niente : 
M, 1.131 2.236 
Pr’ 0 l 
5 0 0 0.5 l 0 
E qi 1.5 0.4 0.6 1.0 0.4 1.1 
I ( 9.5 5.8 18 3.6 5.8 ee 





maximum error of the corrected velocity profile, per 
cent of (1%; — us) 

Er maz = Maximum error of the uncorrected Taylor profile, 

per cent of (u, — u» 
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which is Taylor’s weak-shock solution. The higher-order equa- 
tions for uv” are linear and of the form 


du /dx — 2AuYy™ = W,(x) (10) 


where W,,(x) is a complicated function of lower-order terms. The 
second- and third-order coefficients, in which the fluid proper- 
ties, Pr”, 7, and s, appear explicitly are obtained by solving Eq. 
(10) and are 


u’?) = (B — 1) sech? Ax(log cosh Ax) (11) 


6?) = — (y — 1)(1 — B) sech? Ax(log cosh Ax) + (y — 1) X 
[D sech? Ax — 1]/2 (12) 


u = (1+ F,)sech? Ax tanh Ax[(log cosh Ax)? — log cosh Ax] + 
sech? Ax[F,Ax + F; tanh Ax] (13) 


go) = (¥ — 1) (1 + F,) sech? Ax tanh Ax(log cosh Ax)? — 


(y — 1) FeAx sech? Ax + sech? Ax X 
tanh Ax|F, log cosh Ax + F;] (14) 


The constants appearing in Equations (11), (12), (13), and (14) 
are 


B = s(y + 1) + 4A%y — 1)(1 — 1/Pr”)/Pr"(y + 1); 
F, = B(B — 2) 


D = (2 — Pr”)/((Pr” — 1)/y + 1); 


Fy, = (vy — 1)(8 — 1)(B + D — 1) 


F, = s(y — 1)(8B/2 — (2s + 17 — 1)/2 + 


2A4/Pr" — D/2 — 1] — 2yBA/(y — 1+ Pr”) 
F; = s(y — 1)(1/2 — 3B/2 + (2s + 1)(y — 1)/2 + 
(y+ — 1)D/2y — 2yA(D/y + 1y/(y¥ + 1)) + 
3yBA/(y — 1+ Pr”) 
F, = —(y — 1)} Fy + (yvA/Pr")[s(y — 1) — B — (D — 1)]} 


As convergence could not be proved, the corrections e?u\?) + 
eu, and ¢76°?) + «39 were compared to the exact corrections 
as calculated from analytical solutions® in the cases indicated 
in Table 1. 
locity profiles, which exceeds the temperature-profile error, is 


The maximum error in the corrected Taylor ve 


also given. The exact and series corrections for Pr” | ee 
0.5, are shown in Fig. 1. 

The slow convergence of the series for uw” decreases the accuracy 
of the correction when s > 1/2. Interestingly the Taylor solu- 
tion becomes more accurate as the viscosity-temperature ex- 
ponent, s, increases from 0 to 1. 

In the intervalO < Pr” < 1, the constants B, F,, Fs, ete., vary 
over a much greater range of values than in the interval 1 < 
Pr” < « When Pr” = 0, the series appears to diverge when e > 
(y — 1)/24 
of the Navier-Stokes shock structure as un” — O, such behavior 


From Gilbarg’s® discussion of the limit behavior 


is not surprising. 
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Jet-Compressor Efficiencies as Influenced 
by the Nature of the Driving and Driven Gases 


Suzanne S. Eichacker and Harold J. Hoge* 

Quartermaster Research and Engineering Center, U.S. Army, 
Natick, Mass. 

April 2, 1960 


HE EFFICIENCY OF A JET COMPRESSOR (ejector) is highest 

when the driving fluid is a gas of high molecular weight and 
heat capacity, and the driven fluid is a gas of relatively low 
molecular weight and heat capacity. This is shown by Fig. 1 
A previously reported research! on the system air:air has been 
extended to eight other systems. The maximum efficiency 
attained previously was 0.206; the maximum in the present in- 
vestigation is 0.489, reached with the system helium: Freon—113 
This is considerably higher than any previously reported efli- 
ciency for a jet compressor. The CO,:CO2 curve and the air: air 
curve (the latter reproduced from the previous report) show that, 
with a given ejector, the efficiency of self-entrainment does not 
change much from gas to gas. 

Each point plotted in Fig. 1 represents the highest efficiency 
attainable at a fixed entrainment ratio (m/MM/), when the outlet 
pressure is varied. The results of 71 separate runs are represented 
in the figure. Efficiencies were computed from the formula 


n = (m/M)(ho — hz)/(hy — hw), 


wherein 7 is efficiency, m and M are the respective mass rates of 
flow of the driven and driving fluids, hy-h; is the isentropic en- 
thalpy rise of unit mass of driven fluid, and h;—/h is the isen- 
tropic enthalpy drop of unit mass of driving fluid. The enthalpy 
changes are computed between the actual initial states and the 
observed final (outlet) pressure. The calculated efficiencies de 
pend only on initial and final states and are independent of any as 
sumptions that may be made about the mixing process. The fact 
that the jet compressor delivers the two fluids in a mixed con 
dition is ignored in this definition. 

The notation, experimental procedures, and methods of re 
ducing data used in the present note are the same as those used 
in the previous report. This includes the methods of measuring 
temperatures, pressures, and flow rates, and the method of locat- 
ing efficiency maxima. One gas (Freon—113) could not be passed 
through a dry gas meter, and so was metered by the driving-fluid 
nozzle. Nozzle position remained as in the previous experiments 
Driving-fluid velocities at the nozzle throat and exit were calcu 
lated as in the previous report. The ratio of exit velocity to 

* The authors wish to acknowledge the valuable assistance of Dr. David 
L. Fiske in various phases of the work now being reported 
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Fic. 1. Maximum efficiency versus entrainment ratio for various 
fluid combinations. The driven fluid is given first. 
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Fic. 2. Maximum efficiency versus molecular-weight ratio. 


throat velocity (/;*) had the following values: helium, 1.39; 
air, 1.43; CO,, 1.48; Freon—12, 1.46; Freon—113, 1.52 

Fig. 2 shows the dependence of maximum efficiency on the 
molecular-weight ratio of the fluids employed (driven/driving). 
Each point plotted in Fig. 2 is taken from the peak of the corre- 
sponding curve in Fig. 1. The ordinate is designated n»2 to indi- 
cate that the plotted efficiencies are the highest that could be 
obtained when two variables (outlet pressure and entrainment 
ratio) were independently varied. The equation given in the 
figure was derived to represent the data. The numerical values 
of the constants have significance only for the particular appara- 
tus used. The point of highest efficiency (He: Freon—113) shows 
the greatest departure from the straight line. This is probably 
to be taken as a warning that molecular-weight ratio alone is 
insufficient to determine mixing behavior. One would expect 
the isentropic exponent k Cp/C, to be involved. The observa- 
tion that deviates most widely from the line involves helium, 
which has the highest value of k (1.667) of all the gases employed; 
and Freon—113, which has the lowest value (1.078). The relation 
between efficiency and molecular-weight ratio shown by Fig. 2 
is of considerable practical importance. It shows that steam, 
the driving fluid most commonly used in jet compressors, can 
be expected to give relatively low efficiencies because of its rather 
low molecular weight. 

It seemed plausible that efficiency might be related to the 
velocities of the two streams as they entered the mixing channel. 
It was possible to make a rough calculation of these velocities for 
five of the points plotted in Fig. 2. Driving-fluid velocities were 
obtained as described above; driven-fluid velocities were com- 
puted from the observed pressure drops, assuming isentropic 
expansion from the initial stagnation temperature. When 72 was 
plotted versus the ratio of the velocities of the two streams, a 
strong dependence was found. As velocity-ratio (driven/driving ) 
increased from 0.21 to 0.34, efficiency increased from 0.177 to 
0.359. The dependence appeared to ~be linear over the short 
range covered by the data. Stagnation temperatures of both 
driving and driven fluids were in all cases room temperature, 
which was taken to be 5380°R. Ona few occasions measurements 
were made with the thermocouples provided, to see if there was 
any heat gain or loss to the apparatus. Heat transfer was found 
to be small enough to be neglected 

The goals of high efficiency and high compression ratio desired 
in practice are not completely compatible. It was found that 
the highest efficiencies were always associated with low com- 
pression ratios. The highest compression ratio corresponding 
to any of the points plotted in Fig. 1 was 2.67, and most points 
had compression ratios below 1.5. This may be due to the fact 
that in our apparatus the ratio of the mixing-tube channel area 
to the driving-fluid nozzle-exit area was large (6.97). Work is 
under way to determine the effect of reducing this ratio to 


about 3. 
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An Approximate Equation for the ‘‘Choke- 
Line’’ of a Compressor 


G. T. Csanady 


The University of New South Wales, 
School of Mechanical Engineering, Kensington, Australia 


March 28, 1960 


Manrses and Lewis? have pointed out the similarity be 
tween the pressure-ratio-versus-inlet-mass-flow-coefficient 
characteristic of a steam or gas turbine and the analogous char- 
acteristic of an expansion (Laval) nozzle. An extension of the 
idea to a compressor and a compression nozzle leads to an ap- 
proximate expression for the ‘‘choke line’ of the compressor. 
as shown below 
It is easily shown from the equations governing the flow in 
Laval nozzles that if the throat velocity is sonic in a compression 
nozzle (i.e., if the nozzle is choked) and if the flow is isentropic 
(in which case the nozzle is simply a reversed expansion nozzle), 
exhausting to a reservoir in which the density is p, speed of sound 
a», the outlet mass-flow coefficient is constant with varying ps and 
dy as 


G/p2a2 = const (1 


where G = mass flow. 

The ‘“‘choke line”’ in the performance map of a compressor con- 
nects operating points at which the flow is nearly isentropic, 
the choke just setting in. It is, therefore, not surprising to 
discover that Eq. (1) approximately holds on the choke line 
To show the correspondence for an actual case, consider the 
results of Carter, ef a/.' Outlet density or temperature have not 
been published, and therefore, a polytropic change will be as- 
sumed between inlet and outlet, with an exponent » = 1.5 (small 
stage efficiency = 85 per cent). Then 
+1)/2n » 


p22 = pidi(p2/ pi)" 


as is easily shown by elementary methods. The equation of the 


choke line becomes 


G=K Pi A;( po/ pr )0-853 (3 
with K = const. This represents a line with slight upward 
curvature, leaving the abscissa at a certain G; = A p; a), when 


po/pr = Fitting the constant to the data, the comparison with 
the measured data, Eq. (1) becomes 
p2/ 10 1.47 2.33 3.13 3.4 4.25 
choking flow Ib 

sec. (observed ) 22 32 49 60.5 65 75 
calculated 

(fit at p2/pr: = 

1.47) 23 32 47 60 62.5 77.5 


The agreement is clearly satisfactory, and could possibly be 
improved by considering variations in efficiency in the recalcu- 
lation according to Eq. (2). Many other actual choke lines have 
a shape roughly in accord with Eq. (3) 
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Determination of Pitch and Yaw Rates 
of Models in Free-Flight Ranges 
by Magnetic Induction 


Bernard S. Golosman* 


National Aeronautics and Space Administration, 
Ames Research Center, Moffett Field, Calif. 


April 7, 1960 


Yr CARDS and shadowgraph pictures are presently used to 
determine pitch and yaw of ballistic models in free-flight 
ranges. Both methods yield information at only a limited num- 
ber of points. The magnetic induction methodt here outlined 
yields continuous information on pitch and yaw rates, from 
which continuous records of pitch and yaw angles can be obtained. 
The apparatus for determining rate of pitch consists of two 
long, rectangular multiturn loops, or coils, of wire. One coil is 
located just above and the other just below the test range with 
the length and width of the loops corresponding to the length and 
width, respectively, of the test range (see Fig. 1). The model, 
containing a permanent bar magnet oriented parallel to its 
axis, is fired from a gun down the test range between the upper 
and lower coils. The magnet generates a voltage in the coils 
proportional to the rate of change of pitch angle of the model. 
This voltage is amplified and plotted against time, giving a 
record of pitch rate of the model at any point along its path of 
flight. A similar pair of coils along the sides of the test range 
measures the yaw rate. An experimental setup of the system 
was made in an existing free-flight range at the Ames Research 
Center. 

If the model is far enough from the ends of the coil so that end 
effects can be neglected, if the model is moving along the center 
line of the coil system, and if the pitching and yawing rotations 
are small in amplitude, then 


E, = K,(d6/dt) 
and 
E, = K.(dy/dt) 


where @ and y are the pitch and yaw angles, respectively, E,, is the 
voltage induced in the pitch coil as a result of pitch motion, d6/dt, 
and E, is the voltage induced in the yaw coil. Thus from a 
record of Ep, versus time, approximate values of @ may be ob- 
tained by integration (given an initial value of 6 from a shadow 
graph station). Similarly, one may obtain the approximate 
values of y. 

More precise equations for EF, and E, include an additional 
factor which is a function of @and y. (See reference 1.) There 
are several secondary effects neglected in the above equations 
which may introduce errors, such as (1) change in induced 
voltage due to the model being off the range center line, (2) 
cross-coupling between pitch and vaw channels when the model is 
off center line, (3) coil end effects, and (4) voltages due to com- 

* Aeronautical Research Scientist. 

+ This method was originally suggested by Mr. John Dimeff of Ames 


Research Center. 





Fic. 1. Schematic layout of experimental setup. 
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Fic. 2. Comparison of coil and shadowgraph data for angle of 
attack (= angle of pitch). 


bination of model roll and magnet misalignment with the model 
axis. 

These errors have been analyzed (reference 1) and found to 
be small enough that by proper attention to the design of the 
system, or by correction for the errors, they need not limit the 
accuracy and applicability of the method. 

More troublesome sources of errors, less amenable to analysis, 
are the noise voltages due to extraneous magnetic fields. The 
apparatus had to be designed to minimize the effect of these fields 
To block 60-cycle noise voltages, the coils were mounted in a long 
mild steel circular pipe and a bucking voltage was applied. The 
remaining harmonics of 60 cycles (5 per cent of the signal ampli 
tude) were recorded just before the actual run and subtracted 
from the signal record on the assumption that they were con 
stant during the time of the run. 

Voltages are induced in a coil when it moves in a magnetic 
field. Coil motion can be caused by gun recoil, muzzle blast, 
sabot impact, and bow waves from the model. The effect of 
bow waves from the model was reduced by a felt-lined cardboard 
The effects 
of gun recoil, muzzle blast, and sabot impact were reduced by the 


tube between the model path and the coil assembly 


use of shock mounts and blast shields. 

The strongest magnetic field was due to the earth. The steel 
pipe used to reduce the 60-cycle noise reduced the earth’s mag- 
netic field by two to one. Under these conditions the noise 
voltage due to coil vibration was 7 per cent of the signal voltage 
This vibration noise was the chief limitation on accuracy. 

By aligning the range with the horizontal component of the 
earth’s magnetic field and cancelling the vertical component 
with Helmholtz coils, one should be able to reduce further the 
vibration noise to 0.2 per cent of the signal; however, it was not 
practical to go to these lengths in this case. To avoid inter 
ference with the Helmholtz coils, the 60-cycle shielding pipe 
would need to be of nonmagnetic material 

Conical models were fired down the range and oscilloscope 
records were obtained of the pitch and yaw coil voltages. The 
voltage-versus-time records were converted to records of pitch 
and yaw angles by integration, using values of K predetermined 
by rotation of the model in the coil prior to the run. In the 
absence in this case of a shadowgraph station at the beginning 
of the coil, initiai values for the angles of pitch and yaw were not 
available. Since for conical models the variations of pitch and 
yaw angle were known to be sinusoidal, their derivatives should 
have been maximum when the angles were equal to zero. There 
fore, the voltage peaks in the records were assumed to correspond 
to zero angles of pitch and yaw. 

The method described in reference 2 and values of drag co 
efficient from reference 3 were used to convert the records of 
angles of pitch and yaw versus time to records of angles versus 
distance. Each record was then approximated by a sine wave 
The sine waves were projected past Station No. 1 for comparison 
with the shadowgraph data as shown in Fig. 2. The correlation 
between the projected pickup coil data and the shadowgraph 
data was fairly good. 


REFERENCES 
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The Effect of Thrust and Drag Load 
on the Aeroelastic Behavior of Booster Systems 


John E. Stevens 
Chief of Space Mechanics, Astronautics Division, 
Chance Vought Aircraft, Incorporated 


April 4, 1960 


INTRODUCTION 
Lo slender, multistage booster systems are subjected to 
severe aeroelastic problems while they are traversing the 
lower atmosphere. These problems are associated (1) with the 
joad magnification which results from lateral structural deforma- 
tions that occur due to the angle of attack which is imposed on 
the vehicle due to maneuvers, wind shear, and gusts, (2) with 
control problems which occur due to the changes in stability 
derivatives due to aeroelastic effects, and (3) with the basic 
booster-system, control-systems stability problem. This par- 
ticular development is concerned with the effective reduction in 
lateral stiffness of the booster system which results from the 
column-type loadings on the vehicle. Although a well designed 
booster system will not fail as a column under the action of drag 
and thrust, the work done by these forces tends to reduce the 
apparent lateral stiffness of the booster system and hence in- 
crease the aeroelastic load magnification on a vehicle subjected 
to angle of attack and reduce the vibration frequencies of the 
booster system. This problem is likely to be more severe on 
solid-rocket multistage boosters with high acceleration and hence 
large design dynamic pressures and high acceleration forces than 
it is on large liquid booster systems with lower accelerations and 
lower design dynamic pressures 
EQUATION DEVELOPMENT 
The equations defining the elastic behavior of a missile or 
booster system can be written making use of Lagrange’s equa- 
tions. In most eases, for the lateral deformation of a booster 
system, it is convenient to write an inertia matrix [17] and an 
influence coefficient matrix [C] which define the inertial and 
structural characteristics of a particular vehicle. If we use ma- 
trix notation it is then possible to write the kinetic energy and 
potential energy in a structure as follows: 
2K.E. = {q}'[el'/Wllelig (1) 
and 
2P.E. = {q}'lel'IC]“lelias (2) 
where the lateral displacement of the structure is defined as a 
linear combination of deflection shapes which include a lateral 
translation mode and a rigid-body pitch mode as well as elastic 
modes. Naturally the coefficients of the stiffness matrix for the 
rigid-body modes will be zero. The lateral displacement of the 
panel points along the booster system is then defined by the ma- 
trix equation 
16} = [elias (3) 
where {6} is a column matrix, [g] is a rectangular matrix with as 
many rows as there are panel points to define the lateral displace- 
ment of the booster and as many columns as there are generalized 
shapes selected to define the booster rigid-body displacement and 
deformation shapes, and {g} is a column matrix with as many 
elements as there are generalized coordinates. Associated with 
Eq. (3) we can write a matrix product which defines the slope of 
the deformed structure in the in-flight direction at each panel- 


point location. 
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}a} = [6}} 4} (4) 

We may then write the generalized forces due to the transverse 

pressure distribution, if we assume that the lift per unit length 
along the missile can be written as 

dL/dx = f{(1/2)pV?, M, qi, 2, Qs, Gir Ge, 


The virtual work done by the lift can be written as (the vehicle 


(5) 


has a length /) 
AW = {4q}{Q} = tad’ f, { o(x)}(dL/dx) dx 


Hence the generalized forces are 


F 
{Q} = f } o(x)} (dL /dx) dx (6) 


The equations which define the behavior of the vehicle when 
thrust and drag effects are neglected then become 


: 6a) = ! dL 
[eo] MW) [ela + [el’(C]“lelig} = | o(x)} dx (7) 
0 dx 


The solution of Eq. (7) can be simplified by using orthogonal 
modes as generalized shapes. 

If we now wish to consider the effects of the thrust and drag 
forces, we consider the work done on ihe vehicle due to these 
external forces. At any point along the vehicle at a particular 
time the compressive force, F,(x), along the vehicle’s longitudinal 


axis can be written as (x = 0 at the nose of the vehicle) 
l *dCp r —_ (1 2)pV? S Cp . dm 
FAzx) = - : pV?S - dx + dx 
2 go & M go ax 
where 
Cp = total drag coefficient 
(1/2)oV2 = dynamic pressure 
Ss = a reference area 
7 = the thrust 
M = the total mass of the vehicle 
(dC p/dx) = the change in the drag coefficient per unit length along the 
vehicle 
(dm/dx) = the mass per unit length along the vehicle 


We can now write the virtual work done by the thrust and drag 


forces as 
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aWrv = taal’ (ff, (Gob dx + Ther tort’) tal 


and 
‘ ! > = = ognf = 
1Or.v} = (f, F.(x)|6\'|6| dx + Ti ¢r\\6r}’) ig} (8) 


where, (O}fg} is defined in Eq. (4), T is the thrust from the rocket 
engine, |@,}’{q} is the slope of the deflection mode of the booster 
at the rocket nozzle, and {g7}’{q} is the lateral displacement of 
the booster at the rocket nozzle. 

We can now introduce Eq. (8) into Lagrange’s equation along 
with Eqs. (1), (2), and (6), and define the equations for the 
behavior of the booster system with the column effects included. 


/ 
fel’(Wleliaf + (terici '[y] a F.(x){6 |’ (6| dx — 
0 


l 
; dL 
rer\ ir} ) ig) = | e(x)} x (9) 
0 dx 


For the static aeroelastic problem we are usually satisfied with 
the particular solution to the differential equations in Eq. (9). 
When these values of the generalized coordinate |g} are obtained, 
it is possible by making use of Eq. (5) and the inertia distribution 
to calculate the shear and bending moment in the missile shell 
and the flexible-airframe aerodynamic derivatives for the par- 
ticular flight condition being investigated. 

DISCUSSION 

The method defined in this paper has been used for the design 
of a particular four-stage solid-propellant booster system. The 
effects of the column loading on this particular vehicle were sig- 
nificant. In a particular flight condition the change in bending- 
moment distribution resulting from including column-loading 
effects is shown in Fig. 1. Although the effects on loadings were 
important, the effects of column loading on flexible-body aero- 
dynamic-center location were of even more importance. It was 
concluded from this study that the column-loading effects should 
not be neglected in the design of booster systems which are long 
and slender, particularly when the booster system was subjected 


toa high dynamic pressure. 


Empirical Equation for the Wake-Spreading 
Limits Behind Three-Dimensional Bodies 


W. J. Billerbeck, Jr. 
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[Is PRACTICAL aerodynamic design problems of various types, 
it is often of value to calculate the approximate wake limits 
behind bodies of various shapes. The equations for wake spread- 
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Fic. 1. Variation of wake width with drag coefficient for sharp- 
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ing behind two-dimensional bodies have been discussed by 
Prandtl! and others, and an equation for the wake limits with the 
constants evaluated for airfoils of various types in subsonic flow 
However, 


is presented (reference 2)—viz., ¢ = 0.68 Ca'/*e' 
much less information appears to be available on wakes behind 
three-dimensional bodies. A plot of a very limited amount of 
data available to the writer on bluff bodies is presented here. 
The data on the disc and the square plate were obtained from 
reference 3, while the data on cones were obtained from unpub- 
lished tests at Reynolds numbers of 10° to 10°. It appears that 
an equation of the same form as that developed for the two- 


dimensional bodies roughly fits the data shown. This equation is 


¢ = 0.68 Cy'/*e/ 


Where ¢ is the wake half-width in terms of body diameters, 
Ca is the drag coefficient based on the maximum cross-sectional 
area and e is the distance aft of the trailing edge of the body in 
body diameters. 

It seems particularly interesting, and perhaps significant, that 
the constant derived for three-dimensional bodies is identical to 
that obtained for the wake spreading behind two-dimensional 
bodies. It would appear worth while to conduct some addi- 
tional experiments to provide a more complete confirmation of the 
equation at low values of Ca, and to include data on configu- 
rations which are not sharp-edged for purposes of comparison. 
It may be discovered that this equation will also apply at super- 
critical Reynolds numbers to bodies which are not sharp-edged. 
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simplified by the author’s provision of many subheadings. 

Matrer UsuaLLty DELETED: Photographs or illustrations of 
little technical interest and not showing advances in general 
practice. Too detailed tabular matter (general results of such 
tables may be included in the text). Lengthy descriptions of 
materials or processes or of preliminary experiments or theories 
that preceded final results; salient features only are of interest. 

REFERENCES: References should be numbered consecutively 
and grouped together at the end of the manuscript, with only 
the corresponding number being mentioned in the text. The 
arrangement should be as follows. For books: 'Durand, W. F., 
Aerodynamic Theory, 1st Ed., Vol. 1, p. 23; Julius Springer, 
Berlin, 1934. For magazines: ‘England, C. R., Crawford, 
A. B., and Mumford, W. W., Some Results of a Study of Ultra- 
Short-Wave Transmission Phenomenon, Proc. IRE, Vol. 20, 
No. 12, pp. 481-482, March, 1933. Please give author, title, 
edition, volume, number, page, publisher, and place and date of 
publication as indicated. Omission of one required fact causes 
much extra editorial work and possible inaccuracies. 


ILLUSTRATIONS: Illustrations should accompany manuscripts, 
and each should always be referred to in the text by number. 
Drawings or graphs should not be larger than 12 X 16 inches, 
and must be made with jet black India ink on white paper or 
tracing cloth, the latter being preferred. Do mot use typewriter 
for lettering. The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but it should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
should be drawn in with black ink, omitting the finer divisions. 
In the case of finely ruled paper, only blue-lined paper can be 
accepted. Tracing paper and blue y Sm are not acceptable. 
Lettering and all markings must be large enough to be readable 
after reduction to single-column width (35/\,in.). Mail rolled or 
flat; never fold. Drawings that cannot be reproduced (including 
pencil drawings) will be returned to the author for redrawing, 
thus delaying publication of the paper. Photographs should be 
distinct and show clear black and white contrasts. They must 
be on glossy white paper. Avoid round or oval photographs. 


CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph they should be placed below and well outside 
the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use ‘‘Fig. 1” (not Figure 1), “Figs. 3 and 4,” etc., in both the text 
and the numbering of illustrations. In the text, “Eq. (1)” or 
“Eqs. (1) and (2)” is used, not “Equation (1).” In captions, 
legends, and in table headings, write all words in full; do not 
abbreviate, except for “Fig.”” and “‘Eq.” 


MATHEMATICAL WoRK: Formulas may be typewritten or 
carefully written in pen and ink, the writing to be large enough so 
that it can be marked for the printer. Considerable space for mark- 
ing should be allowed above and below allequations. All compli- 
cated equations should be repeated on separate sheets with plenty 
of space left for marking. The solidus should be used for sim- 
ple fractions appearing within the text. Make all expressions 
clear to the typesetter. Greek letters used in formulas should 
be clearly designated by name in the margin of the manu- 
script. The difference between capital and lower-case letters 
should be clearly distinguished and care taken to avoid confu- 
sion between zero (0) and the letter (0), between the numeral 
(one) and the letter (ell) and the prime (’), between alpha and 
a, kappa and k, u and mu, v and nu, n and eta. All sub- 
scripts and exponents should be clearly distinguished. Avoid 
complicated exponents and subscripts. Dots and bars over 
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When it is necessary to repeat a complicated expression, it should 
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in the American Standard Association ‘‘Letter Symbols 
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carefully checked Standard abbreviations should be used, and 
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